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Abstract 


We  present  a  semi-analytical,  semi-numerical  method  to  calculate  the  diffraction  of 
elastic  waves  by  an  irregular  topography  of  arbitrary  shape.  The  method  is  a  straightforward 
extension  to  three-dimensions  of  the  approach  originally  developed  to  study  the  diffraction  of 
SH  waves  (Bouchon,  1985)  and  P-SV  waves  (Gaffet  and  Bouchon,  1989)  by  two-dimensional 
topographies.  It  relies  on  a  boundary  integral  equation  scheme  formulated  in  the  frequency 
domain  where  the  Green  functions  are  evaluated  by  the  discrete  wavenumber  method.  The 
principle  of  the  method  is  simple.  The  diffracted  wavefield  is  represented  as  the  integral  over 
the  topographic  surface  of  an  unknown  source  density  function  times  the  medium  Green 
functions.  The  Green  functions  are  expressed  as  integrals  over  the  horizontal  wavenumbers. 
The  introduction  of  a  spatial  periodicity  of  the  topography  combined  with  the  discretization 
of  the  surface  at  equal  intervals  results  in  a  discretization  of  the  wavenumber  integrals  and  in 
a  periodicity  in  the  horizontal  wavenumber  space.  As  a  result,  the  Green  functions  are 
expressed  as  finite  sums  of  analytical  terms.  The  writing  of  the  boundary  conditions  of  free 
stress  at  the  surface  yields  a  linear  system  of  equations  where  the  unknowns  are  the  source 
density  functions  representing  the  diffracted  wavefield.  Finally,  this  system  is  solved 
iteratively  using  the  conjugate  gradient  approach.  We  use  this  method  to  investigate  the 
effect  of  a  hill  on  the  ground  motion  produced  during  an  earthquake.  The  hill  considered  is 
120m  high  and  has  an  elliptical  base  and  ratios  of  height-to-half-width  of  0.2  and  0.4, along  its 
major  and  minor  axes.  The  results  obtained  show  that  amplification  occurs  at  and  near  the 
top  of  the  hill  oyer  a  broad  range  of  frequencies. ,  For  incident  shear  waves  polarized  along 
the  short  dimension  of  the  hill  the  amplification  at  the  top  reaches  100%  around  lOHz  and 
stays  above  50%  for  frequencies  between  1.5Hz  and  20Hz.  For  incident  shear  waves 
polarized  along  the  direction  of  elongation  of  the  topography,  the  maximum  amplification 
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occurs  between  2Hz  and  5Hz  with  values  ranging  from  50%  to  75%.  The  results  also  show 
that  the  geometry  of  the  topography  exerts  a  very  strong  directivity  on  the  wavefield 
diffracted  away  from  the  hill  and  that  at  some  distance  from  the  hill  this  diffracted  wavefield 
consists  mostly  of  Rayleigh  waves. 
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Introduction 


Studies  of  the  effect  of  the  irregular  topography  of  the  earth  surface  on  seismic  waves 
have  been  motivated  by  two  types  of  observations.  The  first  investigations  (Tatel,  1954; 
Gilbert  and  Knopoff,  1960;  Key,  1967;  Hudson,  1967;  Greenfield,  1971)  were  stimulated  by 
observations  of  anomalous  arrivals  on  seismograms  recorded  in  or  near  regions  of  rough 
topography.  These  investigations  led  to  the  identification  of  these  arrivals  as  energy  scattered 
by  the  topographic  features  and  consisting  mostly  of  Rayleigh  waves.  Subsequent 
investigations  (  Boore  et  al.,  1981;  Bannister  et  al.,  1990;  Gupta  et  al.,  1990;  Wagner  and 
Langston,  1992;  Clouser  and  Langston,  1995)  have  confirmed  these  findings  and  have  shown 
that  topographic  scattering  is  an  important  source  of  seismic  coda. 

Another  observation  that  lent  importance  to  the  problem  was  the  recording  of  extremely 
high  accelerations  (at  the  time  the  highest  ground  accelerations  ever  recorded)  at  a  site 
located  on  a  topographic  ridge  during  the  San.  Fernando,  California,  earthquake  of  1971. 
Numerical  simulations  carried  out  for  the  particular  configuration  of  the  site  concluded  that 
the  accelerations  recorded  had  been  amplified  by  the  local  topography  (Boore,  1972; 
Bouchon,  1973).  The  interest  originally  incited  by  the  San  Fernando  earthquake  recordings 
has  since  continued.  Field  experiments  have  been  performed  on  hills  and  mountain  ridges 
(Davis  and  West,  1973;  Rogers  et  al.,  1974;  Griffiths  and  Bollinger,  1979,  Tucker  et  al.,  1984; 
Umeda  et  al.,  1986;  Carver  et  al.,  1990;  Pedersen  et  al.,  1994a;  Nechtschein  et  al.,  1995)  and 
they  have  confirmed  that  amplification  -  sometimes  quite  large  -  occurs  near  the  top  of  a  hill 
or  mountain  and  along  the  crest  of  a  ridge. 

Some  visual  observations  made  in  epicentral  areas  of  strongly  felt  earthquakes  such  as  the 
presence  of  churned  and  shattered  earth  at  the  top  of  prominences  (Muir,  1912;  Hadley,  1964; 
Nason,  1971)  and  the  disruptions  of  rocks  and  boulders  near  hill  crests  (Oldham,  1899;  Clark, 
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1972;  Bolt  and  Hansen.  1977;  Uraeda,  1990)  also  indicate  the  occurrence  of  intense  shaking 
in  elevated  areas  of  rugged  topography.  Further  interest  in  this  problem  was  generated  by 
observations  of  the  damage  patterns  of  recent  earthquakes  like  the  1987  Whittier  Narrows, 
California,  earthquake  (Kawase  and  Aki,  1990)  and  the  1989  Loma  Prieta  earthquake  (Ponti 
and  Wells,  1991;  Hartzell  et  al.,  1994)  and  by  recordings  made  during  the  1994  Northridge, 
California,  earthquake  of  extremely  high  accelerations  (1.78g,  one  of  the  highest  acceleration 
ever  recorded  during  an  earthquake)  at  a  site  located  on  top  of  a  hill  (Shakai  et  al.,  1994). 

Motivated  by  all  these  observations  and  by  the  need  to  include  and  quantify  these  effects 
in  seismic  risk  assessment  and  microzonation  studies,  a  remarkable  amount  of  work  has  been 
done  on  the  theoretical  investigation  and  on  the  numerical  simulation  of  the  diffraction  of 
seismic  waves  by  irregular  topographies  (Trifunac,  1973;  Sabina  and  Willis,  1975;  Wong  and 
Jennings,  1975;  Singh  and  Sabina,  1977;  Sills,  1978;  Sanchez-Sesma  and  Rosenblueth,  1979; 
England  et  al.,  1980;  Ilan  and  Bond,  1981;  Bard,  1982;  Sanchez-Sesma  et  al.,  1982;  Shah  et 
al.,  1982;  Wong,  1982;  Ohtsuki  and  Harumi,  1983;  Sanchez-Sesma,  1983;  Zahradnik  and 
Urban,  1984;  Bard  and  Tucker,  1985;  Bouchon,  1985;  Geli  et  al.,  1988;  Kawase,  1988;  Gaffet 
and  Bouchon,  1989;  Sanchez-Sesma  and  Campillo,  1991;  Pedersen  et  al.,  1994b;  Takemiya 
and  Fujiwara,  1994).  Laboratory  model  studies  of  topographic  scattering  observed  in  the  field 
have  also  been  conducted  (Rogers  et  al.,  1974;  Anooshehpoor  and  Brune,  1989). 

The  theoretical  and  numerical  results  obtained  for  ridge  type  topographies  consistently 
predict  the  occurrence  of  amplification  at  the  ridge  crest.  The  levels  of  amplification 
calculated,  however,  almost  systematically  underestimate  the  actual  amplifications  observed 
in  the  field.  Several  hypotheses  have  been  advanced  to  account  for  this  discrepancy  (Griffiths 
and  Bollinger,  1979;  Bard,  1982;  Bard  and  Tucker,  1985;  Geli  et  al.,  1988;  Hartzell  et  al., 
1994).  Among  them  is  the  two-dimensionality  of  the  topographic  geometries  assumed  in  the 
calculations.  Recent  field  experiments  using  seismic  arrays  and  conducted  at  the  site  where 
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high  accelerations  were  recorded  during  the  Northridge  earthquake  emphasize  the  three- 
dimensional  nature  of  the  seismic  response  of  the  hill  (Lee  et  al.,  1994;  Spudich  et  al.,  1995). 
The  only  three-dimensional  theoretical  investigation  reported  in  the  geophysical  litterature  to 
date  seems  to  be  the  study  by  Sanchez-Sesma  (1983)  of  scattering  by  an  axi-symmetric 
surface  irregularity. 

We  present  in  this  paper  a  relatively  simple  semi-analytical  semi-numerical  approach  to 
calculate  the  diffraction  of  elastic  waves  by  an  irregular  three-dimensional  topography  of 
arbitrary  shape.  The  method  relies  on  a  boundary  integral  equation  formulation  and  on  the 
use  of  the  discrete  wavenumber  method  to  evaluate  the  Green  functions.  It  constitutes  the 
straightforward  extension  to  the  3D  case  of  the  method  originally  developed  to  study  the 
diffraction  of  SH  waves  (Bouchon,  1985)  and  of  P-SV  waves  (Gaffet  and  Bouchon,  1989)  by 
2-dimensional  topography. 

Description  of  the  method 

We  consider  the  problem  of  a  seismic  wavefield  impinging  from  below  on  the  surface  of 
an  elastic  homogeneous  half-space.  We  denote  by  a,  P  and  p  the  P  and  S  wave  velocities  and 
density  of  the  medium  and  by  X  and  p  its  Lame  coefficients.  We  define  a  cartesian  coordinate 
system  (x,y,z)  in  which  z  is  chosen  to  be  the  downward  vertical  direction.  Then,  the 
displacement  field  m  at  a  point  x  may  be  written  in  the  form: 

M,(  X  )  =  (  X  )  -f- 1  Qj{s)Gij{  X  ,\s)ds  (1) 

s 

where  u°  is  the  incident  source  field,  Qj{s)  is  an  unknown  source  density  function,  Gij  is  the 
medium  Green  function,  S  denotes  the  diffracting  surface,  and  the  subscript  i  indicates  the 
component  of  displacement  considered. 

Each  point  s  of  the  surface  acts  like  a  diffracting  source  of  strength  Qj{s)  that  may  be 
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represented  as  consisting  of  3  forces  acting  along  the  (x.y,z)  directions  (Coutant,  1989). 

In  order  to  perform  the  integration  expressed  by  equation  (1)  we  introduce  a  spatial 
periodicity  in  the  x  and  y  directions  and  discretize  the  surface  at  equal  Ax  and  Ay  intervals. 
We  then  first  assume  that  the  incident  wavefield  has  a  time  dependence  and  express  the 
Green  functions  Gy  in  the  form  of  horizontal  wavenumber  integrals  (Lamb,  1904).  The 
spatial  periodicity  results  in  a  discretization  of  the  wavenumber  integrals  (Bouchon  and  Aki, 
1977;  Bouchon,  1979)  while  the  spatial  discretization  at  constant  Ax  and  Ay  intervals  implies 
a  periodicity  in  the  horizontal  wavenumber  domains  (Bouchon,  1985;  Campillo  and  Bouchon, 
1985).  Equation  (1)  then  becomes: 

IV,  N, 

£/,(  X )  =  t/f  ( X )  +  2  E  Qjih^y)Gij{  X ,  Xy,/^)  Ax  Ay  (2) 

i>i  (,,=1 

with,  using  the  displacement  potential  expressions  given  in  Bouchon  (1979), 
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and  where  Lx=NxAx  and  Ly-NyLy  denote  the  spatial  periodicity  lengths  in  the  x  and  y 
directions,  and  Nx=lMx+\  ,  N^=2My+\ ,  implying  that  andiVy  are  odd  integers. 

The  corresponding  stress  wavefield  may  be  written 

N,  N, 

<5pq{  X  )  =  C!°q{  X  )  +  X  E  QAk^y)Hpqj{  X  ,  X;^,/^)  ^X  Ay  (4) 

4=1  4=1 

with;  ^pqj  ~ 

where  5p^=l  for  ,  and  5p,=0  for  pj^q  ,  and  where  Gpj^q  denotes  the  derivative  of  Gpj 
with  respect  to  the  q  coordinate. 

The  unknown  source  densities  Qj(ix,iy)  are  obtained  by  solving  the  stress-free  boundary 
conditions  at  each  discretized  point  (JxJy)  of  the  surface: 

'^pUxyjy)  ~  ^qUxJy)  JxJy^  ~ 


8 


where  n  denotes  the  unit  normal  to  the  surface. 


Equation  (6)  is  a  system  of  3N.,Ny  linear  equations  which  may  be  rewritten: 


-Vx 


or 


Qji^X’^y) 


^qUxJy) 
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(8) 


with 

CpjUxJyUxJy)  =  riqUxJy)  ,x  i^jJAxAy  (9a) 

^ppUxJy'JxJy)  —  (9b) 

The  expression  for  Cpp{jx,jy\jxJy)  (equation  (9b))  which  is  independent  of  the  local  normal 
direction  is  derived  by  taking  n.{JxJy)=^  in  equation  (9a)  as  discussed  in  Bouchon  and 
Coutant  (1994). 

The  right  hand  side  of  equation  (8)  represents  (with  a  change  of  sign)  the  incident  stress  at 
the  surface.  Its  evaluation  is  straightforward  and  can  be  obtained  using  the  discrete 
wavenumber  method  (Bouchon,  1979)  if  the  source  is  a  point  source  or  an  extended  source. 
In  the  case  where  the  source  is  a  vertically  propagating  S  wave  of  unit  spectral  amplitude  and 
polarized  along  the  x  direction,  as  will  be  considered  in  the  next  section,  the  expressions  for 
the  incident  displacement  and  stress  wavefields  are: 

,  ■  •  ,  i  i'  ■  >',•'/  '  '.•I  * 

.  0) 

t^°x{x,y,z)=e  P 

.  CD  ^ 

C7^?(j^,y,z)=0  except  a",(^,y,z)=a^(j:,y,z)=zp-^e'p‘'  (10) 
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r"=o 

In  the  case  where  the  source  is  an  obliquely  propagating  plane  wave,  equation  (3  a)  must  be 
replaced  by: 

iKm-r  iKtUy 

If  =ifO  , _ t. .  h  -1,0  , - L  . 

Lx  Ly 

where  k°  and  ky  are  the  horizontal  wavenumber  components  of  the  incident  wave. 

The  solution  of  the  system  of  equations  (8)  yields  the  strength  and  phase  of  the  diffracting 
sources  distribution  Qj(ix,iy)  from  which  the  wavefield  can  be  evaluated  anywhere  in  the 
medium  using  equation  (2).  The  system  is  solved  for  each  frequency  and  the  solution  is 
multiplied  by  the  spectrum  of  the  source  time  function  of  the  incident  wavefield  and  Fourier 
transformed  to  the  time  domain.  The  effect  of  the  periodicity  of  the  structure  is  removed  by 
using  complex  frequencies  as  discussed  in  Bouchon  and  Aki  (1977)  and  Bouchon  (1979). 

The  discretization  intervals  zkx  and  Ay  are  determined  for  each  frequency:  at  low 
frequencies  a  minimum  number  of  discretization  points  Nx  and  Ny  is  assumed,  insuring  that 
the  topographic  surface  is  well  defined;  as  the  frequency  increases  these  numbers  are  defined 
such  that  there  are  at  least  2.5  points  per  seismic  wavelength,  that  is: 


As  a  result,  at  high  frequency,  the  size  of  the  system  increases  like  the  square  of  the 
frequency.  In  order  to  reduce  this  size  we  remove  from  the  matrix  equation  (8)  all  the  terms 
which  lie  below  a  certain  amplitude  threshold  as  discussed  in  the  two-dimensional  case  by 
Bouchon  et  al.  (1995)  and  Schultz  et  al.  (1995),  and  we  invert  the  resulting  sparse  system 
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iteratively  using  the  conjugate  gradient  method.  A  test  showing  the  accuracy  of  the  solution 
obtained  with  this  procedure  is  presented  in  the  next  section. 

Effect  of  a  hill 

We  now  use  the  method  to  investigate  how  the  topography  of  a  hill  affects  the  seismic 
ground  motion.  The  topography  considered  is  displayed  in  Figure  1.  The  base  of  the  hill  is 
an  ellipse  of  minor  and  major  axes  equal  to  600m  and  1200m.  The  top  of  the  hill  stands  at 
120m  and  its  flanks  have  a  cosine  shape.  The  size  of  the  model  considered  (defined  by  the 
periodicity  lengths  in  the  two  horizontal  directions)  is  2000m  by  2000m.  The  medium  P- 
wave  velocity  is  3000m/s  and  the  Poisson  ratio  is  0.25.  We  first  consider  the  case  of  a 
vertically  incident  shear  wave  polarized  along  the  minor  axis  of  the  ellipse.  We  assume  that 
this  incident  wave  has  a  Ricker-pulse  time  dependence  centered  about  lOHz.  The  resulting 
ground  motion  produced  along  the  profile  that  crosses  the  hill  along  its  minor  axis  is  shown  in 
Figure  2.  One  observes  a  clear  amplification  of  the  ground  motion  at  and  near  the  top  of  the 
hill.  The  hill  also  diffracts  part  of  the  incoming  vertically  propagating  S-wave  into  a 
horizontally  propagating  P-wave  and  a  Rayleigh  wave.  Near  the  bottom  of  the  hill  these  two 
diffracted  waves  have  comparable  amplitude,  but,  as  distapce  increases,  the  Rayleigh  wave 
becomes  predominant  as  the  surface  P-wave  attenuates  more  rapidly.  The  ground  motion 
produced  along  the  perpendicular  profile  that  runs  along  the  major  axis  of  the  hill  is  displayed 
in  Figure  3.  The  results  show  that  there  is  almost  no  diffracted  wavefield  radiated  in  this 
direction. 

The  parameters  used  in  this  simulation  are  as  follows:  31  equally  spaced  frequencies 
between  OHz  and  30  Hz  are  considered;  below  15Hz  all  the  matrix  terms  are  included,  while 
above  15Hz  an  amplitude  threshold  of  1%  is  applied.  The  validity  of  the  solution  obtained  by 
removing  the  matrix  terms  that  lie  below  this  amplitude  threshold  may  be  estimated  by 
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comparing  it  to  the  complete  full  matrix  solution  (Figure  4).  This  comparison  is  made  at 
15Hz  for  the  horizontal  ground  displacement  calculated  along  a  profile  running  from  the  top 
of  the  hill  to  the  edge  of  the  model  and  following  the  minor  axis  of  the  topography.  At  this 
frequency  the  sparse  system  retains  only  7%  of  the  non-zero  terms  of  the  full  system. 

Snapshots  of  the  horizontal  ground  displacement  are  presented  in  Figure  5.  They  show 
the  strong  directivity  that  the  shape  of  the  topography  exerts  on  the  diffracted  wavefield;  most 
of  the  diffracted  energy  is  radiated  in  the  direction  perpendicular  to  the  direction  of 
elongation  of  the  topography  and,  as  seen  also  in  Figure  3,  only  insignificant  energy  is 
diffracted  in  directions  close  to  the  direction  of  elongation.  Most  of  the  wavefield  diffracted 
away  from  the  hill  has  been  generated  in  the  summit  area  of  the  hill,  as  is  also  clearly  seen  in 
Figure  2.  There  is  also  some  energy  diffracted  by  the  bottom  of  the  hill  but  it  is  much 
weaker. 

In  order  to  assess  the  amplification  of  ground  motion  produced  by  the  hill  and  investigate 
its  frequency  dependence,  we  computed  the  spatial  variation  of  the  peak  horizontal 
displacement  produced  for  4  frequencies  of  excitation  (Figure  6).  In  each  case  the  source  is  a 
vertically  incident  S-wave  pulse  polarized  along  the  minor  axis  of  the  hill.  In  the  first  case 
(X,5=/i)  the  center  frequency  of  the  pulse  is  such  that  its  wavelength  is  equal  to  the  height 
h  of  the  hill.  In  the  other  cases  the  center  wavelengths  of  the  incident  pulse  are  twice,  three 
times,  and  four  times  the  height.  The  peak  horizontal  displacement  is  inferred  by  combining 
the  time  traces  of  the  two  horizontal  displacement  components.  The  results  show  that 
amplification  occurs  at  and  near  the  top  of  the  hill  at  all  frequencies  considered.  This 
amplification  is  maximum  at  the  top  of  the  hill.  The  zone  of  high  amplification  is  elongated 
and  follows  the  ridge  of  the  elongated  topography.  The  width  of  this  zone  is  about 
proportional  to  the  incident  seismic  wavelength.  The  amplification  reaches  100%  for  Xs=h, 
90%  for  Xs=2h,  70%  for  X,5=3/i,  and  60%  for  Xs=4h.  Some  attenuation  takes  place  on  the 
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flanks  of  the  hill. 


We  now  investigate  the  case  where  the  incident  S-wave  is  polarized  along  the  direction  of 
elongation  of  the  topography.  The  horizontal  ground  displacement  produced  along  the  two 
perpendicular  profiles  crossing  the  hill  is  displayed  in  Figure  7.  Some  amplification  still 
occurs  in  the  summit  area  of  the  hill.  This  amplification,  however,  is  not  as  strong  as  the  one 
produced  when  the  incident  wave  is  polarized  along  the  short  dimension  of  the  topography. 
The  shape  of  the  topography  still  exerts  the  same  directivity  on  the  wavefield  diffracted  away 
from  the  hill,  no  significant  energy  being  radiated  in  the  direction  of  elongation  of  the 
topography.  The  absence  of  Rayleigh  waves  on  the  profile  which  runs  along  the  minor  axis 
of  the  hill  is  due  to  a  node  of  radiation  of  P-waves  in  this  direction. 

Maps  of  corresponding  peak  horizontal  displacements  as  a  function  of  incident 
wavelength  are  presented  in  Figure  8.  Amplification  still  takes  place  at  or  near  the  top  of  the 
hill  for  the  three  longest  wavelengths  considered.  The  maximum  amplification  (75%)  now 
occurs  at  the  longest  wavelength  (ks=4h)  and  diminishes  as  the  frequency  increases  (65%  at 
A.5=3;z,  35%  at  Xs=2h,  0%  at  X.5=/i). 

In  order  to  investigate  more  fully  the  frequency  dependence  of  the  topographic 
amplification,  we  present  in  Figure  9  the  evolution  of  the  displacement  amplitude  ratio 
between  the  top  and  the  base  of  the  hill  as  a  function  of  frequency.  For  incident  waves 
polarized  along  the  elongated  direction  of  the  topography,  the  peak  amplification  occurs 
between  2Hz  and  5Hz.  la  this  frequency  band  the  ground  motion  at  the  top  of  the  hill  is 
amplified  by  values  ranging  from  50%  to  75%.  In  contrast  almost  no  amplification  occurs  for 
frequencies  higher  than  lOHz.  For  motion  transverse  to  the  direction  of  elongation  of  the  hill, 
-the  peak  amplification  occurs  around  lOHz  where  it  reaches  about  100%.  The  amplification 
of  the  transverse  ground  motion  stays  above  50%  over  a  broad  frequency  range  extending 
from  1.5Hz  to  20Hz. 
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Conclusion 


We  have  presented  a  semi-analytical  semi-numerical  method  to  calculate  the  diffraction 
of  seismic  waves  by  a  three-dimensional  topography  of  arbitrary  shape.  The  formulation 
relies  on  a  boundary  integral  equation  scheme  where  the  Green  function  are  evaluated  by  the 
discrete  wavenumber  method.  The  boundary  conditions  of  free  stress  at  the  surface  lead  to  a 
linear  system  of  equations,  which  is  solved  iteratively  by  the  conjugate  gradient  method.  We 
have  used  this  formulation  to  investigate  the  effect  of  a  hill  on  earthquake  ground  motion. 
We  found  that  for  the  configuration  studied  -  the  hill  is  120m  high  and  its  base  is  an  ellipse 
1200m  in  length  and  600m  in  width  -  amplification  occurs  at  and  near  the  top  of  the  hill  over 
a  broad  range  of  frequencies.  Amplification  levels  are  the  highest  for  shear  waves  polarized 
transversely  to  the  direction  of  elongation  of  the  hill.  The  amplification  of  transverse  motions 
at  the  top  of  the  hill  reaches  100%  at  lOHz  and  stays  above  50%  between  1.5Hz  and  20Hz. 
For  motion  along  the  direction  of  elongation  of  the  topography,  the  maximum  amplification 
occurs  between  2Hz  and  5Hz  with  values  ranging  from  50%  to  75%.  At  some  distance  from 
the  hill  the  diffracted  wavefield  consists  mostly  of  Rayleigh  waves.  The  geometry  of  the  hill 
can  exert  a  drastic  directivity  on  the  diffracted  wavefield. 
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Figure  legends 


Figure  1 

Shape  of  the  topography  considered.  The  distances  and  heights  are  expressed  in  meters 
and  the  elevation  contour  lines  are  spaced  at  10m  interval.  The  hill  has  an  elongated  elliptical 
shape  with  an  aspect  ratio  of  two. 

Figure  2 

Horizontal  (top)  and  vertical  (bottom)  displacements  produced  by  a  vertically  propagating 
S  wave  incident  on  the  topography  of  Figure  1.  The  shear  wave  is  polarized  along  the  short 
dimension  of  the  hill  and  the  profile  runs  along  the  same  direction.  A  cross-section  of  the 
topography  is  displayed  on  the  left  of  the  figure. 

Figure  3 

Horizontal  displacement  produced  by  a  vertically  propagating  S  wave  incident  on  the 
topography  of  Figure  1 .  The  shear  wave  is  polarized  along  the  short  dimension  of  the  hill 
while  the  profile  runs  along  the  major  axis  of  the  hill.  The  corresponding  topographic  cross- 
section  is  displayed  on  the  left  of  the  figure. 

Figure  4 

Comparison  in  amplitude  and  phase  between  the  full  and  the  sparse  matrix  solutions  for 
the  horizontal  displacement  at  a  frequency  of  15Hz.  The  profile  considered  runs  from  the  top 
of  the  hill  to  the  edge  of  the  model  and  follows  the  minor  axis  of  the  topography. 

Figure  5 

Snapshots  of  the  horizontal  displacement  produced  for  a  vertically  incident  S  pulse 
polarized  along  the  short  dimension  of  the  topography.  The  component  of  displacement 
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displayed  is  the  one  parallel  to  the  incident  polarization  direction.  The  time  sequence  is 
indicated  and  the  time  interval  between  successive  frames  is  l/32s.  The  incident  pulse  is  a 
Ricker  wavelet  with  a  center  period  of  0.138s,  corresponding  to  a  wavelength  equal  to  twice 
the  hill  height.  The  first  seven  snapshots  are  dominated  by  the  primary  incident  pulse  which 
illuminates  the  topography.  The  elliptic  base  of  the  hill  is  visible  on  snapshot  6.  The  last 
nine  snapshots  show  the  diffracted  wavefield  being  generated  mostly  in  the  summital  area  of 
the  hill  and  propagating  away  from  the  hill. 

Figure  6 

Maps  showing  the  spatial  distribution  of  horizontal  displacement  amplitude  for  four 
different  frequencies  of  excitation.  Xs  denotes  the  shear  wavelength,  and  h  is  the  height  of 
the  topography.  The  dotted  line  indicates  the  base  of  the  hill. 

Figure  7 

Horizontal  ground  motion  produced  along  two  profiles  running  parallel  to  the  major  (top) 
and  minor  (bottom)  axes  of  the  hill  for  a  vertically  incident  S  wave  polarized  along  the 
elongated  dimension  of  the  topography.  The  corresponding  topographic  cross-sections  are 
displayed  on  the  left  of  the  figure. 

Figure  8 

Maps  showing  the  spatial  distribution  of  horizontal  displacement  amplitude  for  four 
different  frequencies  of  excitation,  %  denotes  the  shear  wavelength,  and  ft  is  the  height  of 
the  topography.  The  dotted  line  indicates  the  base  of  the  hill. 

Figure  9  • 

Amplitude  ratio  between  the  ground  motions  at  the  top  and  at  the  base  of  the  hill  in  the 
case  of  incident  shear  waves  polarized  along  directions  parallel  and  transverse  to  the  direction 
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of  elongation  of  the  topography.  The  reference  point  at  the  base  of  the  hill  is  chosen  at  some 
distance  from  the  hill  so  as  not  to  be  affected  by  the  presence  of  the  hill. 
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Amplitude  of  horizontal  displacement 
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Incident  S  wave  polarized  along  the  topography  minor  axis 


Amplitude  of  horizontal  displacement 
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Incident  S  wave  polarized  along  the  topography  major  axis 
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Abstract 


We  present  a  boundary  integral  equation  /  conjugate  gradient  formulation  to  study  the 
propagation  of  seismic  waves  through  complex  geological  structures.  The  method  is  aimed  at 
extending  the  range  of  applications  of  boundary  integral  equations  or  boundary  elements 
methods  to  geological  models  of  relatively  large  size  or  complexity.  We  show  that  the 
system  of  equations  which  expresses  the  boundary  conditions  at  the  medium  interfaces  and 
which  is  inherent  to  the  biem  or  bem  approach  can  be  drastically  reduced  in  size  and  that  only 
10  to  20%  of  the  terms  of  this  system  contribute  significantly  to  the  solution.  The  boundary 
conditions  may  thus  be  expressed  in  the  form  of  a  very  sparse  linear  system  which  can  be 
inverted  iteratively  by  the  conjugate  gradient  method.  We  use  this  approach  to  investigate  the 
effect  of  a  sedimentary  layer  of  varying  thickness  on  local  and  regional  seismic  wave 
propagation.  We  show  that  the  amplitude  of  ground  motion  and  the  duration  of  shaking  are 
drastically  increased  when  the  sediment/bedrock  interface  is  rough  and  irregular  relatively  to 
the  case  where  it  is  flat. 
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Introduction 


In  recent  years  several  methods  using  boundary  integral  equations  or  boundary  elements 
have  been  developed  to  calculate  the  propagation  of  seismic  waves  through  complex 
geological  structures  (e.g.  Sanchez-Sesma  and  Esquivel,  1979;  Dravinski,  1983;  Bouchon, 
1985;  Campillo  and  Bouchon,  1985;  Campillo,  1987;  Bravo  et  al.,  1988;  Kawase,  1988; 
Coutant,  1989;  Kawase  and  Aki,  1989;  Gaffet  and  Bouchon,  1989,  1991;  Papageorgiou  and 
Kim,  1991;  Sanchez-Sesma  and  Campillo,  1991;  Mossessian  and  Dravinski,  1992;  Schultz 
and  Toksoz,  1993;  Pedersen  et  al.,  1994).  One  common  feature  of  these  techniques  is  that 
they  require  the  solution  of  a  linear  system  of  equations  which  expresses  the  boundary 
conditions  and  which,  when  applied  to  real  cases,  can  be  very  large.  The  size  of  this  system 
is  indeed  the  most  limiting  factor  in  the  applicability  of  the  method. 

The  aim  of  the  present  study  is  to  investigate  ways  of  reducing  the  size  of  the  system  of 
linear  equations  inherent  to  the  boundary  integral  equations  formulation.  We  shall  show  that 
this  can  be  achieved  simply  by  removing  from  the  equations  the  Green  function  terms  which 
lie  below  a  cenain  amplitude  threshold  and  which  do  not  contribute  significantly  to  the 
solution.  We  shall  limit  the  scope  of  this  study  to  the  case  of  SH  waves  and  we  shall  present 
an  application  of  the  method  to  regional  seismic  wave  propagation. 

Irregular  surface  topography 

We  begin  the  presentation  with  a  rather  simple  case  which  involves  a  plane  SH  wave 
vertically  incident  on  an  irregular  topography.  Denoting  by  x  and  z  the  horizontal  and 
vertical  axes,  the  reflected  and  diffracted  displacement  wavefield,  at  a  particular  frequency  co, 
may  be  written  in  the  form  of  an  integral  over  the  surface  5  of  a  source  density  function  q 
times  the  medium  Green  function  G  ,  that  is: 
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V{x,z)=jqis)G{x,z-^s^Zs)ds  (1) 
s 

Using  the  discrete  wavenumber  method  to  evaluate  the  Green  function  (Bouchon  and  Aki, 
1977)  and  discretizing  the  surface  at  equal  Ax  interval  along  the  x-axis.  equation  (1)  may  be 
expressed  in  the  form  (Bouchon,  1985): 

n«)=^i;ey  z  ^ ^ -  (2) 

j=l  m=-N/l  * 

with: 


Ym  =  (y  Im  (y„)<0  . 


and  where  p  is  the  shear  modulus,  P  is  the  shear  wave  velocity,  L  is  the  required  periodicity 
length  of  the  surface  topography,  and  where  the  Qj  are  the  unknown  strengths  of  the  antiplane 
forces  located  at  the  discretized  points  iXj,Zj)  of  the  surface.  These  forces  represent  the 
sources  of  the  scattered  wavefield.  In  this  formulation,  the  surface  is  assumed  to  have  a 
periodic  shape.  In  practice  this  periodicity  length  (like  the  size  of  the  model  in  finite- 
differences  methods)  is  chosen  large  enough  so  that  repeated  topographies  do  not  affect  the 
time  domain  solution.  N  is  the  number  of  discretized  points  and  is  required  to  be  an  odd 

integer.  The  discretization  interval  Ax  =  -^  is  set  according  to  the  wavelength  considered. 

The  unknovim  strengths  Qj  of  the  forces  are  determined  by  matching  the  boundary 
conditions  of  free  stress  at  each  one  of  the  surface  points.  In  order  to  write  down  the 
corresponding  equations,  let  us  denote  by  V/y  and  (T,y  the  displacement  and  surface  stress 
produced  at  a  point  i  of  the  surface  by  the  surface  force  distribution  element  Qj  acting  at  point 
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j.  We  have,  using  the  discrete  wavenumber  formalism  of  equation  (2)  and  denoting  by  and 
n,  the  X  and  z  components  of  the  normal  to  the  surface: 


Qj  N/Z 


S 

m=-N/l 


dVij  dVn 

(4, 


that  is: 


'm=-N/2 


with  sgnizi-Zj)=l  forz,>zy  ,  =-l  for  z,<zy . 


When  i=j  the  derivation  expressed  by  equation  (4)  is  undefined  and  equation  (5a)  must  be 


replaced  by  (Bouchon  and  Coutant,  1994): 


Let  us  now  consider  the  case  where  a  vertically  monochromatic  plane  SH  wave  is  incident 
upon  the  surface.  The  incident  displacement  and  stress  at  a  point  i  of  the  surface  may  be 
written  in  the  form: 


.  CO 

y  •  =  e  p 
^01  ^ 
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The  boundary  conditions  require-  that  the  stress  associated  with  the  scattered  wave  field 
cancells  the  incident  stress  along  the  surface,  that  is: 

N 

XcJy  =-ao,- ,  for/=l,iV  (7) 

M 

which  may  be  written  in  matrix  form: 


CQ  =  F 


(8) 


with: 


K 

Ym 


for  iitj 


Cii  = 


LK 

2  L 


Fi  =  -Ooi  • 


The  amplitude  of  the  C  matrix  elements  is  largest  for  the  diagonal  terms  and,  on  average, 
decreases  with  increasing  [i—j ]  separation. 

In  order  to  investigate  whether  all  the  terms  of  the  matrix  are  necessary  to  yield  an 
accurate  solution,  we  performed  the  inversion  of  the  system  after  retaining  only  the  high 
amplitude  terms  of  the  matrix.  This  makes  C  a  sparse  matrix  and  mostly  retains  the  terms 
closest  to  the  diagonal.  We  then  use  the  conjugate  gradient  method  to  invert  the  now  sparse 
linear  system.  A  description  of  the  conjugate  gradient  method  and  the  corresponding 
software  can  be  found  in  Press  et  al.  (1994).  A  comparison  of  the  results  obtained  with  the 
corresponding  full  matrix  inversion  is  presented  in  Figure  1 .  The  topography  considered  is  a 
200m  high  hill  with  a  cosine  shape  and  a  base  length  of  800m.  The  shear  wave  velocity  is 
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lOOOm/s.  A  vertically  propagating  Ricker-type  SH  pulse  with  5Hz  center  frequency  is 
incident  on  the  surface.  The  periodicity  length  chosen  is  4000m.  The  calculation  is  done  for 
a  time  window  length  T  of  2s  and  over  32  frequencies  starting  at  OHz  and  spaced  at  equal 

A/  =  -^  interval.  Characteristic  of  the  discrete  wavenumber  method  and  later  removed  from 

the  time  domain  solution,  a  constant  imaginary  part  equal  to  -i—  is  added  to  the  angular 

frequency.  The  number  of  discretized  points  A/^  is  equal  to  201. 

The  traces  displayed  in  the  bottom  half  of  Figure  1  are  the  results  of  the  full  matrix 
inversion.  They  correspond  to  receivers  located  along  the  topography,  from  the  top  of  the  hill 
to  2000m  away  from  the  summit.  These  seismograms  are  compared  with  those  shown  in  the 
upper  half  of  the  figure  which  are  obtained  by  retaining  only  15%  of  the  matrix  elements  C^. 

In  this  case  the  removal  of  elements  is  achieved  using  a  simple  amplitude  criterion:  Only  the 
15%  elements  with  the  highest  amplitude  are  kept.  The  15%  threshold  used  here  is  somewhat 
arbitrary  but  it  shows  that  a  good  solution  is  reached  by  retaining  only  the  high  amplitude 
terms  of  the  matrix.  •  Doing  so  reduces  drastically  the  size  of  the  system  of  equations  and 
results  in  a  considerable  gain  in  memory  space  and  computer  time. 

Laterally-varying  layered  medium 

Let  us  now  consider  the  case  of  a  multilayered  medium  where  the  shape  of  the  interfaces 
is  irregular.  If  we  denote  by  Q‘{1)  and  Q^{1)  the  unknown  force  distributions  respectively  at 
the  top  and  at  the  bottom  of  layer  /,  and  by  S\l)  and  S^{1)  the  eventual  source  displacement- 
stress  vectors  incident  respectively  on  the  top  and  bottom  interfaces  of  layer  Z,  the  boundary 
conditions  of  continuity  of  displacement  and  stress  across  the  interface  separating  layer  / 
(above)  from  layer  Z-i-1  (below)  may  be  written: 

Zaf;(/)  Q){1)+  e;^(/)  +  5f(/)=  SaJ(/+l)  gyCZ+D  +  ^RZ+l) 

y=i  j=i  j=i  7=1 
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for  i-\,2N 


(9) 


where  the  aijil)  are  the  displacement  and  stress  Green  functions  for  medium  /.  The 
superscript  bt  indicates  that  the  point  of  evaluation  of  the  field  (first  index)  is  at  the  bottom  of 
the  layer  while  the  point  of  application  of  the  force  (second  index)  is  at  the  top  of  the  layer, 
and  so  on. 


Thus,  using  the  discrete  wavenumber  formalism  described  above,  we  have: 


1  N/l  -iKn(Xi-Xj) 

2;  i - 

m=-N/2 


\  N/l  lc„ 

X  K  ~  +  ^8^  (Zi-Zj)n^.  for  iVy 


i=-M/2  Tm 


fli+w.KO  =  “ytt 


fori=l,A^:;=l,iV 


(10) 


with: 

Im(Yj<0 

sgn{zi-Zj)=l  for  Zi^j,  =-l  for Zi<Zj,  fora“,a*' 
sgn{zi-Zj)=l  for Zi>Zj,  =-l  forz,<ry,  fora*^,a'* 

Let  us  consider  the  case  where  the  free  surface  is  flat  and  located  at  z=0.  The  boundary 
conditions  at  the  bottom  of  the  first  layer  may  then  be  expressed: 

+  Sf{l)  =  Za‘(j{2)Q)i2)+  ^afjC2)Q){2)  +  S\{2)  (11) 

y=i  7=1  7=1 
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where  dij  are  the  half-space  displacement  and  stress  Green  functions  for  medium  1,  that  is,  in 
discrete  wavenumber  form: 


4/(1)  =  ^  I  ^ 


1  N/1  k  -  k 

4+M;(1)  =  Z  [K-  ZT'^sgn  iZi-Zj)n^.  +  (n^  ^ 

m--Nn  ifn  \m 


for/=l,^^,  j=\,N. 


Equations  (9)  and  (11)  lead  to  a  linear  system  of  the  form: 


CQ  =  F  (12) 


where,  denoting  by  nl  the  number  of  layers: 


A**(l)  -A  "(2)  -A '*(2) 


0  A^'(2)  A**(2)  -A“(3)  -A'^(3)  0  0 


0  A%/-1)  A^^nl-l)  -A"(n/) 


5'(2)-^^(l) 
5‘(3)-5^(2)  ,  „ 


The  dimension  of  the  linear  system  expressed  by  equation  (12)  is  2iVx(n/-l).  The  number  of 


discretized  points  N  representing  each  interface  varies  with  the  frequency  considered.  At  low 
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frequencies,  this  number  is  set  to  a  minimum  value,  and,  as  the  frequency  increses,  its  value  is 
chosen  such  that  there  are  at  least  2.5  points  per  seismic  wavelength,  that  is: 


A/'  =  2.5  X  L  X 


_ frequency _ 

minimum  wave  velocity 


(13) 


In  general,  a  different  value  of  N  is  defined  for  each  interface,  according  to  the  minimum 
shear  wave  velocity  in  the  surrounding  layers.  The  2.5  points  per  seismic  wavelength 
criterion  arises  from  tests  of  accuracy  performed  for  the  configurations  considered  in  this 
study. 

Because  of  the  linear  dependency  of  N  on  distance  and  frequency,  numerical  simulations 
of  high  frequency  seismic  wave  propagation  over  long  distance  ranges  cannot  be  attempted 
through  a  straighforward  inversion  of  equation  (12),  as  the  inversion  of  the  coniplete  system 
using  classical  methods  would  be  prohibitive.  Thus  we  proceed  as  for  the  topography  case. 
After  normalization,  we  retain  only  the  largest  terms  of  the  Green  function  matrix  C  and  use 
the  conjugate  gradient  method  to  solve  the  resulting  sparse  system.  The  configuration  studied 
is  displayed  at  the  bottom  of  Figure  2.  The  medium  consists  of  a  three-layers  crust  overlying 
a  mande  half-space.  In  this  simulation  we  investigate  the  effect  that  a  highly  irregular 
subsurface  topography  between  sediments  and  bedrock  has  on  ground  motion.  For  this 
reason  we  assume  that  the  other  two  interfaces  are  flat  although  they  are  similarly  treated  in 
the  boundary  integral  equation  scheme.  The  source  is  a  line  of  shear  dislocation  across  a 
vertical  fault  plane  perpendicular  to  the  laterally-varying  structure  and  located  at  a  depth  of 
10km.  The  source  displacement-stress  wavefield  incident  on  the  top  and  bottom  interfaces  of 
the  source  layer  4,  expressed  in  diiscrete  wavenumber  form,  is  given  by: 


5P(4)  =  V„  = 


m=-N/2  Tm 
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and 


(14) 


,  h  ‘ 


N/2 


r  l^m[nx—+sgniZi-Zo)n,jQ 


^“Ym  ' 


m=-N/2 


where  is  the  source  seismic  moment  and  {Xg,Zo)  denote  the  source  coordinates. 

Equations  (14)  are  obtained  by  derivation  of  equations  (3)  and  (5a)  with  respect  to  x^,  after 
replacing  the  index  y  of  the  diffracting  source  by  the  primary  source  index  o. 

In  order  to  make  sure  that,  like  for  the  topography  case,  the  size  of  the  system  can  be 
reduced  and  still  yield  a  good  solution,  we  compare  in  Figure  2  the  surface  displacement 
obtained  by  full  matix  inversion  with  the  one  calculated  by  retaining  only  the  largest  terms  of 
the  matrix.  We  choose  as  a  criterion  -  rather  arbitrarily  -  to  eliminate  all  the  matrix  elements 
which  have  an  amplitude  smaller  than  1%  of  the  amplitude  of  the  largest  term  in  the  row 
(which  is  the  diagonal  element).  We  perform  this  comparison  at  a  frequency  of  2Hz,  which  is 
the  median  value  of  the  frequency  range  considered. 

As  the  agreement  between  the  two  solutions  is  satisfactory,  we  extend  this  1  %  amplitude 
threshold  criterion  to  all  frequency  calculations  above  2Hz.  Below  2Hz  we  choose  to  retain 
the  full  matrix  as  the  size  of  the  system  is  then  small  enough  to  carry  out  the  full  system 
inversion.  The  results  obtained  are  showed  in  Figure  3a.  They  display  the  horizontal  ground 
velocity  in  the  direction  tangential  to  the  plane  of  propagation.  The  source  has  a  step-like 
time  dependence  and  the  frequency  range  considered  extends  from  OHz  to  4Hz.  The 
periodicity  length  is  L=  275km.  The  number  of  elements  considered  at  each  frequency  is 
showed  in  Figure  4.  At  low  frequencies  a  minimum  number  corresponding  to  a  minimum 
discretization  interval  Ax  of  1km  designed  to  provide  an  accurate  representation  of  the 
interfaces  geometry  is  assumed.  The  abrupt  change  of  curve  beyond  2Hz  corresponds  to  the 
application  of  the  threshold  criterion.  At  the  highest  frequency  only  12%  of  the  non-zero 
elements  of  the  matrix  are  utilized  to  calculate  the  solution. 
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In  order  to  assess  the  effect  of  the  irregular  sediments-basement  topography  on  the  ground 
motion  we  also  simulated  the  case  where  this  interface  is  flat.  The  results  are  displayed  in 
Figure  3b.  The  configuration  is  the  same  as  the  one  simulated  in  Figure  3a  except  that  the 
sediment  thickness  is  constant  and  equal  to  2km.  The  calculation  is  performed  with  the  same 
parameters  as  in  Figure  3a,  but  with  the  plane  layers  reflection  and  transmission  coefficients. 
The  differences  between  the  two  simulations  are  quite  dramatic.  The  change  from  a  flat 
basement  topography  to  an  irregular  one  considerably  enhances  the  ground  motion, 
particularly  the  duration  of  ground  shaking.  The  irregular  shape  of  the  interface  diffracts  the 
organized  incoming  seismic  energy  into  rays  which  travel  in  all  directions.  Many  of  these 
rays  stay  trapped  in  the  sediments  where  they  give  rise  to  surface  waves  which  propagate 
laterally  across  the  sedimentary  structures  as  has  been  shown  in  several  numerical  simulation 
experiments  (Bard  and  Bouchon,  1980a,b;  Hill  and  Levander,1984;  Levander  and  Hill,  1985; 
Benz  and  Smith,  1988;  Chavez-Garcia  and  Bard,  1989;  Hill  et  al.,  1990;  Papageorgiou  and 
Kim,  1991;  Gaffet  and  Bouchon,  1991;  Frankel  and  Vidale,  1992;  Graves  and  Clayton,  1992; 
Kawase  and  Sato,  1992;  Frankel,  1993;  Jongmans  and  Campillo,  1993).  One  such  train  of 
Love  waves  can  be  clearly  seen  propagating  forward  in  the  basin  closest  to  the  source,  then 
being  reflected  at  the  edge  of  the  basin  and  propagating  backward  across  it  before  being 
reflected  again  at  the  other  edge.  Likewise,  much  of  the  surface  seismic  motion  over  the 
basins  can  be  attributed  to  Love  waves.  Between  the  basins,  in  the  regions  of  bedrock 
outcrops,  the  seismic  motion  is  greatly  reduced,  before  the  amplitudes  pick  up  again  in  the 
next  basin.  This  pattern  is  strikingly  similar  to  the  one  observed  by  Hanks  (1975)  and  Liu 
and  Heaton  (1984)  and  subsequently  studied  by  Vidale  and  Helmberger  (1988)  on  recordings 
of  the  San  Fernando  earthquake  made  across  the  San  Fernando  and  the  Los  Angeles  basins, 
although  in  that  case  Rayleigh  waves  play  the  role  held  by  Love  waves  in  our  study. 

The  case  where  the  deeper  crustal  interfaces  are  also  irregular  is  displayed  in  Figure  5. 
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The  velocity  and  density  of  the  layers  as  well  as  the  shape  of  the  sediment/bedrock  interface 
are  the  same  as  in  the  previous  example,  but  the  upper  crustal  discontinuity  and  the  Moho  are 
now  also  irregular.  All  three  interfaces  are  represented  by  the  same  number  of  discretized 
points  as  in  the  previous  case  and  the  same  amplitude  threshold  criterion  is  applied.  The 
resulting  surface  motion  resembles  the  one  calculated  in  Figure  3a,  indicating  that  the  surface 
motion  is  much  more  sensitive  to  the  topography  of  the  sediment/bedrock  interface  than  to 
the  shape  of  the  deeper  crustal  and  crust/mantle  discontinuities.  In  order  to  assess  the  effect 
of  removing  the  low-amplitude  terms  of  the  matrix  system  for  this  configuration,  we  compare 
in  Figure  6  the  seismograms  obtained  by  inverting  the  full  matrix  with  those  obtained  by 
retaining  all  the  matrix  terms  up  to  IHz  and  removing  the  ones  below  a  1%  amplitude 
threshold  between  IHz  and  4Hz.  Here  again  the  agreement  is  quite  satisfactory  for  many 
applications. 


Conclusion 


We  have  tried  to  extend  the  range  of  applications  of  seismic  wave  propagation  problems 
which  can  be  studied  using  boundary  integral  equation  formulation  by  reducing  the  size  of  the 
linear  system  which  expresses  the  boundary  conditions  and  which  is  usually  the  factor  which 
limits  the  extent  and  the  complexity  of  the  geological  models  which  can  be  investigated.  We 
have  shown  that  only  the  10  to  20%  highest  amplitude  terms  of  the  system  are  important  to 
obtain  a  good  solution.  The  removal  of  the  other  terms  does  not  significantly  affect  the 
results.  The  inversion  of  the  very  sparse  system  is  then  performed  iteratively  using  the 
conjugate  gradient  approach.  As  an  application  of  the  method,  we  have  shown  the  drastic 
effect  that  an  irregular  sediment/bedrock  interface  has  on  the  amplitude  and  duration,  of 
seismic  ground  shaking. 
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Figure  legends 


Figure  1 

Ground  motion  produced  by  a  SH  pulse  vertically  incident  on  a  hill-shaped  topography. 
The  upper  half  of  the  figure  corresponds  to  the  sparse  matrix  solution  while  the  lower  half 
displays  the  full  matrix  solution. 

Figure  2 

Source-medium  configuration  considered  and  corresponding  amplitude  and  phase  of  the 
surface  displacement  at  2Hz.  The  source  location  is  indicated  by  the  star.  The  curves 
(amplitude)  and  symbols  (phase)  display  the  comparison  between  the  full  matrix  solution  and 
the  1%  amplitude  threshold  solution. 

Figure  3 

(a)  Ground  horizontal  velocity  produced  by  an  antiplane  dislocation  for  the  crustal  model 
of  Figure  2.  The  epicentral  distance  range  extends  from  0  to  lOOkm.  A  reduced  velocity  of 
3600m/s  is  applied  to  the  traces,  (b)  Same  as  (a)  when  the  sediment/bedrock  interface  is  flat. 

Figure  4 

Number  of  elements  considered  as  a  function  of  frequency  for  the  calculation  of  Figure 

3a.  :  ,  V  . 


Figure  5 

Ground  horizont^  velocity  produced  by  an  antiplane  dislocation  and  corresponding 
crustal  model.  The  epicentral  distance  range  extends  from  0  to  100km.  A  reduced  velocity  of 
3600m/s  is  applied  to  the  traces. 

Figure  6 


53 


Comparison  of  ground  velocity  time  histories  calculated  by  inverting  the  full  matrix  with 
those  obtained  after  removing  the  terms  which  lay  below  a  1%  amplitude  threshold  for 
frequencies  between  IHz  and  4Hz.  The  crustal  structure  geometry  is  the  one  depicted  in 
Figure  5.  Only  the  shape  of  the  upper  crustal  discontinuities  is  shown  here. 
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ABSTRACT 

In  this  paper,  two  different  methods  to  solve  scattering  problems  in  acoustic  or  elastic  media  are 
coupled  to  enhance  their  usefulness.  The  multiple  multipole  (MMP)  expansions  are  used  to  solve 
for  the  scattered  fields  in  homogeneous  regions  which  are  possibly  unbounded.  The  finite  element 
(FE)  method  is  used  to  calculate  the  scattered  fields  in  heterogeneous  but  bounded  scatterers. 
As  the  MMP  method  requires,  the  different  regions  and  methods  are  coupled  together  in  the.  least 
squares  sense.  For  some  examples,  the  scattered  fields  are  calculated  and  compared  to  the  analytical 
solutions.  Finally,  the  seismograms  are  calculated  for  a  scattering  problem  with  several  scatterers, 
and  complex  geometries.  Thus,  the  hybrid  MMP-FEM  technique  is  a  very  general  and  useful  tool 
to  solve  complex,  two-dimensional  scattering  problems. 
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INTRODUCTION 


Wave  scattering  problems  have  been  investigated  by  different  techniques.  Analytical  solutions  to  the 
integral  equations  do  generally  not  exist  except  for  some  very  simple  geometries.  Analytical  mode 
expansion  is  limited  to  geometries  such  as  circular  cylinders  or  spheres  where  the  modes  decouple 
(Pao  and  Mow,  1973).  Therefore,  numerical  schemes  seem  to  be  the  most  direct  procedure  for 
arbitrary  geometries.  Numerical  boundary  integral  techniques  (Schuster  and  Smith,  1985),  the  T- 
matrix  method  (Waterman,  1969,  1976)  and  MMP  expansions  (Hafner,  1990;  Imhof,  1995a,b)  are 
examples  thereof.  Unfortunately,  they  all  depend  on  either  Greens  functions  or  other  solutions  to  the 
wave  equation  which  tend  to  be  hard  or  impossible  to  find  for  heterogeneous  or  anisotropic  media. 
Thus,  these  methods  are  normally  limited  to  scattering  between  homogeneous  scatterers  embedded 
in  a  homogeneous  background.  As  an  advantage,  these  methods  do  not  encounter  problems  with 
unbounded  domains.  No  artificial  radiating  boundary  conditions  have  to  be  enforced.  In  fact,  the 
scattered  fields  can  be  evaluated  anywhere. 

In  cases  where  the  medium  is  heterogeneous,  finite  element  (FE)  (Zienkiewicz,  1977;  Schwarz, 
1988;  Murphy  and  Chin-Bing,  1989,  1991;  Marfurt,  1984)  or  finite  differences  (FD)  (Marfurt,  1984; 
Kelly  et  ai,  1976;  Virieux,  1986)  techniques  are  routinely  used  to  calculate  the  scattered  wavefields. 
Opposed  to  the  boundary  methods  mentioned  priorly,  FE  and  FD  encounter  serious  problems  with 
unbounded  domains.  The  domain  has  to  be  truncated  and  radiating  boundary  conditions  have  to  be 
enforced.  Even  if  the  domains  are  bounded,  they  are  limited  due  to  computer  memory  and  runtime 
considerations.  For  many  problems  the  distance,  between  inhomogeneities,  source  and  receivers  are 
rather  large  and  thus  result  in  prohibitive  computation  times  and  memory  requirements. 

Many  scattering  problems  exist  which  fall  in  between  these  two  classes.  These  problems  involve 
heterogeneous  regions  which  are  bounded  and  embedded  in  a  homogeneous  background.  Therefore, 
there  is  an  obvious  interest  in  combining  methods  for  unbounded,  homogeneous  domains  with 
methods  which  can  handle  heterogeneous  regions  of  limited  extent  (Su,  1983;  Dubus,  1994).  In 
the  present  paper,  such  a  combination  is  made  between  the  MMP  expansions  and  the  FE  method 
(FEM). 

We  will  apply  the  hybrid  technique  to  both  acoustic  and  elastic  in-plane  scattering  problems 
where  one  or  multiple  heterogeneities  are  embedded  in  a  homogeneous  medium.  Both  source  and 
receiver  are  in  the  homogeneous  region.  All  the  ideas  presented  will  also  hold  if  the  source  and 
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receiver  are  located  in  the  heterogeneity.  We  could  then  use  the  combined  MMP-FEM  technique 
to  construct  the  radiating  boundary  condition.  In  this  work,  we  will  not  investigate  this  usage  of 
the  technique.  Also,  we  will  neglect  the  case  of  anti-plane  wave  motion  (SH)  because  it  can  easily 
be  derived  from  the  acoustic  case. 

This  paper  is  structured  as  follows;  First,  we  will  review  both  MMP  expansions  and  the  finite 
element  method  for  the  acoustic  case.  Next,  we  combine  the  methods  for  the  acoustic  case.  Then, 
we  review  MMP  and  FEM  in  the  elastic  case  and  present  the  combination  thereof.  Finally,  we 
discuss  some  details  of  the  implementation  on  a  computer,  present  solutions  to  some  scattering 
problems  and  compare  them  to  analytical  solutions  where  available. 

ACOUSTIC  THEORY 

We  would  like  to  model  how  an  incident  wavefield  of  angular  frequency  oj  scatters  from 

an  object.  The  situation  is  depicted  in  Figure  1.  The  scatterer  is  heterogeneous  and  embedded 
in  a  homogeneous  background  QP  .  For  the  sake  of  clarity,  we  will  suppress  the  time  factor 
in  all  following  expressions.  Where  necessary,  the  superscripts  ^  and  ^  will  denote  quantities 
which  belong  to  the  homogeneous  outside,  lie  on  the  boundary  between  the  domains  or  are  inside 
the  heterogeneous  region,  respectively.  Quantities  marked  with  a  tilde  are  either  transformed 
quantities  (e.g.,  LU  decomposed)  or  local  quantities  for  one  particular  little  element  where  the 
context  allows  to  infer  the  correct  meaning. 

Homogeneous  Regions:  Multiple  Multipole  Expansions 

In  a  homogeneous  region  expansions  for  the  pressure  fields  are  made  with  exact  solutions  to 
the  homogeneous  wave  equation 

jo 

pO(x)  =  ;^p9pO(x)  ,  where  ^  VPf  :  (V^ -|- fc^)pP  =  0  (1) 

.  j=i 

where  ko  =  uj/ao  is  the  wave  number  and  ao  the  wave  velocity  in  the  homogeneous  region.  The 
factors  pP  are  complex  valued  weighting  coefficients  for  the  different  expansion  functions  .  As 
the  name  of  the  method  implies,  several  multipole  solutions  centered  at  different  positions  are  often 
used  as  expansion  functions.  The  reason  to  use  multiple  multipole  expansions  is  their  local  behavior 
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and  thus  their  ability  to  model  wavefields  scattered  from  complex  geometries  (Imhof,  1995a). 


M  N 

=  E  E  -^H  (ko  |X  -  X^l) 

m=l  n=—N 


(2) 


The  function  is  the  Hankel  function  of  the  first  kind  and  order  n  radiating  outward.  Each 
summation  over  the  index  n  builds  up  one  multipole.  To  enhance  the  convergence,  M  different 
expansion  centers  located  at  are  used.  Since  the  Hankel  functions  have  a  singularity  at  their 
origin,  the  centers  of  expansions  Xm  may  not  be  located  in  the  homogeneous  region  For  each 
expansion  center,  all  orders  between  —N  <n<  +N  are  used  as  basis  functions. 

However,  additional  expansion  functions,  such  as  plane  waves  or  other  special  modes,  can  be 
included.  As  a  result,  MMP  expansions  have,  in  general,  a  smaller  number  of  unknowns  than  com¬ 
parable  methods.  Equations  for  the  weighting  coefficients  are  obtained  by  enforcing  boundary 
conditions  for  the  pressure  and  the  normal  displacement  on  discrete  matching  points  mi  on  the 
boundaries  between  domains.  The  boundary  conditions  between  two  domains  and  are 


5]pfpP(mi)-bP‘-(mi)  = 


^ni-£p?VPf(mi)  +  ^n..VP‘-(mi)  =  ^ni  ■  VP^(m.) 


(3) 


(4) 


where  it  is  assumed  that  the  only  source  is  the  incident  field  P*”^  in  the  homogeneous  domain 
Thus,  we  can  build  the  following  linear  equation  system 


IJO  J  \ 

where  the  submatrices  and  Ufj  contain  Pj^  and  p^fii  •  VPj^  evaluated  at  the  matching  points 
mi.  In  general,  expansions  of  the  form  (1)  are  not  orthogonal.  Thus,  more  matching  points  than 
needed  are  used  and  the  resulting  overdetermined  linear  system  (5)  is  solved  in  the  least  square 
sense  minimizing  the  overall  error  in  the  boundary  conditions.  Due  to  their  definition,  the  matrices 
and  U®  are  rectangular  and  dense  matrices. 
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Heterogeneous  Regions:  Finite  Elements 

Neglecting  source  terms,  waves  propagating  in  an  heterogeneous  region  are  governed  by  the 
general  Helmholtz  equation. 

V  •  VP}  +  p-^k^  P  =  0)  (6) 

where  p  =  p(x)  and  k  =  k{x)  denote  density  and  wave  number,  respectively. 

To  solve  this  equation,  we  partition  heterogeneous  domain  into  small  and  nonoverlapping 
elements  (Zienkiewicz,  1977;  Schwarz,  1988;  Murphy  and  Chin-Bing,  1989).  Commonly,  one 
chooses  triangular  or  quadrangular  elements.  In  each  element,  the  pressure  field  P  is  approximated 
.by  an  interpolation  function.  For  a  quadrangular  element,  the  most  simple  interpolation  function 
to  use  is  the  bilinear  one: 

P{x)  =  Co  +  aix  +  a2Z  +  a^xz  (7) 

Instead  of  directly  using  the  coefficients  aj,  the  polynomial  (7)  is  transformed  into  the  sum  of  simple 
shape  functions  Nj{x)  having  local  support  only.  For  example,  in  a  rectangular  element  of  unit 
size,  Nz{x)  =  xz.  This  shape  function  is  visualized  in  Figure  2.  The  other  ones  are  obtained  by 
rotations  of  —180°,  —90°  and  90°. 

4 

■P(x)  =  (8) 

j=i 

In  fact,  the  complex  valued  weighting  coefficients  pj  are  the  pressure  values  at  the  element’s  corners 
Xi.  They  are  also  known  as  node  points.  The  coefficients  pj  are  called  node  variables.  Also,  the 
interpolation  functions  have  to  satisfy  the  following  conditions: 

Njixi)  =  6ij  (9) 

4 

.^Nj{x)  =  1  for  X  €  (10) 

3=1 

Thus,  the  pressure  in  the  Helmholtz  equation  (6)  is  replaced  by  the  interpolation  (8).  Applying 
Galerkin’s  method,  we  multiply  the  resulting  expression  by  the  test  function  Ni{x)  and  integrate 
over  the  element  Cl. 
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where  the  divergence  theorem  is  used  to  transform  the  first  integral.  The  integrals  are  over  the 
surface  Cl  or  around  the  boundary  f  of  the  element  Cl,  respectively.  Evaluating  (11)  for  all  four 
iVj(x)  will  yield  a  set  of  four  equations  for  the  four  unknown  node  variables  pj. 

These  three  integrals  define  the  local  stiffness  matrix  S,  the  local  mass  matrix  M  and  the  local 
force  vector  f. 


~Sij  = 

If  p-^VNi-VNjdA 

(12) 

JJ  n 

Mij  = 

ff  p-’^k^NiNjdA 

(13) 

JJ  n 

(14) 

For  elements  that  are  in  the  interior,  the  boundary  integral  (14)  is  not  zero,  but  its  contributions 
will  exactly  cancel  with  like  terms  coming  from  neighboring  elements.  One  only  need  to  recall  that 
the  term  is  proportional  to  the  normal  displacement.  But  both  the  normal  displacement 

and  the  pressure  are  continuous  across  boundaries.  Therefore,  only  on  the  domain  boundary  the 
line  integral  has  to  be  taken  into  account  since  it  is  not  cancelled  by  another  term. 

If  the  element  Cl  is  adjacent  to  a  rigid  domain,  the  boundary  integral  (14)  will  vanish,  since 
Ni  =  0.  If  the  element  is  adjacent  to  a  void  domain,  the  integral  (14)  also  vanishes  because 
^  =  0.  In  all  other  cases,  the  boundary  integral  (14)  has  to  be  included.  Assuming  that  ^  can 
be  approximated  by  a  function  similar  to  Ni  along  the  boundary,  we  can  replace  (14)  by 

=  Fij  =  j.p-^NiNjdl.  (15) 

If  the  density  p  and  the  wavenumber  k  are  treated  as  constants  within  each  element,  S,  M  and  F 
can  be  evaluated  exactly.  Once  the  contribution  of  the  various  elements  is  determined,  the  global 
system  of  equations  is  formed  by  mapping  the  local  node  numbers  onto  the  global  node  numbers, 
giving  rise,  to  the  global  pressure  vector  p,  and  combining  all  of  the  subsystems  S,  M  and  f  into 
their  global  counterparts  S,  M  and  f  (Schwarz,  1988).  Both  matrices  S  and  M  are  sparse,  banded 
and  symmetric.  Each  row  of  the  global  matrix  system  can  then  be  reduced  to 

=  0  (16) 
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where  J  is  the  total  number  of  node  variables  or  in  simpler  matrix  form 

K  •  p  -  f  =  0  (17) 

where  K  =  S  —  M  and  the  vector  p  contains  all  the  unknown,  global  nodal  values  pj. 

Coupling  the  Regions 

The  vector  p  obtained  from  the  finite  elements  containing  the  node  variables  can  be  split  into  two 
subvectors  p^  and  p^.  The  node  variables  from  inside  the  domain  are  collected  in  the  vector 

p^.  The  subvector  p^  accommodates  the  node  variables  whose  nodal  points  Xj  lie  on  the  boundary 

dO,^.  Since  the  boundary  dCl^  belongs  to  both  domains,  the  respective  wavefields  and 
have  to  match  along  the  boundary. 

Thus  we  replace  the  node  variables  pf  by 

jo 

pf  =  ^p^P,(xi)  +  P'”‘=(xi)  (18) 

fc=l 

Furthermore,  we  can  also  find  by  evaluating 

^  p^niVPfc(xi)  +  niVP*”‘=(xi)  (19) 


Combining  (16),  (18)  and  (19)  yields  the  hybrid  matrix  system 


where  is  the  total  number  of  node  variables  inside  the  heterogeneous  region  is  the  node 

number  of  the  first  nodal  point  lying  on  the  boundary  The  value  of  is  J^  + 1.  As  before,  J 
is  the  total  number  of  node  points.  Finally,  is  the  total  number  of  functions  used  for  the  MMP 
expansion  of  the  outside  field.  The  complete  system  can  be  written  in  a  more  compact  form  as 


AOI  j  ypO  j  fO  j 


(21) 
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where  is  a  sparse,  diagonally  dominant  and  symmetric  matrix.  Both  and  A^^  are  sparse 
and  rectangular,  while  A®®  is  rectangular,  but  dense  matrix.  The  force  vector  is  sparse,  while 

is  completely  filled. 

The  matrix  A^^  and  the  solution  vector  stem  from  the  interior  problem  which  is  solved  by 
FE.  It  contains  as  many  equations  as  unknowns  and  can  be  solved  exactly.  The  matrices  A^^  and 
A^^  couple  the  interior  problem  to  the  exterior  problem  and  vice  versa.  The  matrix  A^^  and  the 
solution  vector  stem  from  the  exterior  problem  which  is  solved  by  MMP  expansions.  It  has 
to  be  solved  in  the  least  squares  sense  because  there  are  more  equations  needed  than  unknowns 
given.  Therefore,  the  complete  matrix  system  (21)  can  not  be  solved  exactly.  Contrarily,  it  does 
neither  mathematically  nor  physically  make  sense  to  solve  the  complete  system  (21)  in  the  leaist 
squares  sense.  Prior  experience  with  MMP  methods  shows  that  at  least  twice  as  many  equations 
as  unknowns  are  needed  to  obtain  a  reasonable  solution  (Imhof,  1995b).  Unfortunately,  there  are 
in  general  more  unknowns  in  the  interior  than  in  the  exterior.  Thus,  it  is  nearly  impossible  to 
obtain  more  than  twice  as  many  equations  as  unknowns.  Furthermore,  the  solution  in  the  interior 
is  already  an  approximation  to  the  wave  equation.  Solving  the  complete  system  in  the  least  squares 
sense  distributes  the  errors  evenly  over  all  unknowns  which  corrupts  the  solution  in  the  interior 
further. 

Therefore,  the  system  (21)  is  solved  in  two  steps:  first,  the  interior  node  variables  p^  are 
eliminated  by  a  partial  Gaussian  elimination.  Because  the  corresponding  submatrix  A^^  is  derived 
with  the  finite  elements  method  and  thus  diagonally  dominant,  the  Gaussian  elimination  can  be 
performed  without  additional  pivoting. 


The  submatrix  A^^  is  now  an  upper  triangular  matrix.  The  remaining  system  can  then  be  solved 
in  the  least  squares  sense  using  the  normal  equations 


where  the  superscript  ^  denotes  the  complex  conjugate  transpose.  If  desired,  the  values  of  the 
node  variables  p^  are  found  by  back-substitution. 

.  p^  =  .  pO  (24) 
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Practically,  the  system  (21)  is  solved  by  a  combined,  row-wise  LU-QR  algorithm.  Prom  each  new 
row,  the  interior  node  variables  are  Gaussian  eliminated.  Then,  Givens  row  updating  (schwarz, 
1989)  is  performed  on  the  remaining  row.  The  scheme  is  equal  to  normal  Givens  updating  with 
the  first  Givens  rotations  replaced  by  Gaussian  eliminations  instead.  Thus,  the  first  rows  are 
only  LU  decomposed.  All  other  rows  are  additionally  Givens  rotated. 

Remark:  An  Alternative  Solver  Scheme 

Alternatively,  the  system  (21)  can  be  solved  by  the  iterative  scheme: 

A"-p^  =  (25a) 

AOO.pO  ^  fO_AO/.p/_^  (25b) 

Optimally,  each  of  (25a)  and  (25b)  is  also  solved  by  an  iterative  scheme  such  as  the  conjugate 
gradient  method  (Hestenes  and  Stiefel,  1952).  The  first  part  (25a)  is  square,  symmetric  and  sparse. 
The  second  part  (25b)  is  rectangular  and  dense,  but  relatively  small  compared  to  (25a).  The  scheme 
(25)  offers  an  alternative  to  (21),  but  this  has  not  yet  been  tried. 

ELASTIC  THEORY 


Homogeneous  Regions:  Multiple  Multipole  Expansions 

In  a  homogeneous  region  expansions  for  the  displacement  fields  w(x)  =  (u(x),v(x))  are  made 
with  exact  solutions  to  the  homogeneous  wave  equation  (Imhof,  1995b). 


j=l 

(26) 

where 

■  '  wf(x)  -  ^  ■V$,(x)^  ' 

(27a) 

w/(x)  =  Vx{y^j{x)) 

•  (27b) 

and  each  expansion  function  or  satisfies  a  Helmholtz  equation 

iv^  +  kl)^j  =  o 

(28a) 

(v2  -t-  =  0 

(28b) 
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where  ko  and  Iq  are  the  the  wave  numbers  of  the  P-,  respective  S-wave  in  the  homogeneous  region. 
The  factors  <f)j  and  are  complex  valued  weighting  coefficients  for  the  different  expansion  functions 
and  ^j.  Similar  to  the  acoustic  MMP  expansions  (2),  we  choose  multipole  expansions  for 
and 

M  N 

=  E  E  (29a) 

m=l  n=—N 

M  N  ‘  ■ 

=  E  E  (29b) 

m— 1  n=—N 


Equations  for  the  weighting  coefficients  and  ipj  are  obtained  by  enforcing  boundary  con¬ 
ditions  along  discrete  matching  points  m,  on  the  boundaries  between  domains.  The  boundary 
conditions  between  two  domains  and  are  continuity  of  displacement  and  stresses  in  normal 
and  tangential  directions: 


jo 

Ei 

i=i 

jo 

E 

j=l 


,wj(mi)| -|-■«’'“^(In^)  =  u^{mi) 

(30) 

J^;J(mi)| +t;‘"^(mi)  =  v'^(mi) 

CO 

yfii  +  hi  <T”*‘^(nii)  hi  =  hi  <T^(mi)  hi 

(32) 

i|  hi-l-ficr'"‘^(mi)hi  =  fj  cr^(mi)  Hj 

(33) 

mi)  denote  the  stresses  evaluated  at  mf 

due  to  the 

jO 

E« 

j=i 

jO 

E> 

j=i 


displacements  w|’(mi),  w|'(mi),  wj"‘=(mi)  and  wj^(mi),  respectively.  Accordingly,  we  can  build  a 
linear  equation  system 


/  -  U‘"‘^  \ 

_  ^inc 

U  j 

sf  sf  J 

( sf  -  sr  y 

where  the  submatrices  contain  equations  (30-33)  evaluated  at  all  matching  points  nii. 


(34) 
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Heterogeneous  Regions:  Finite  Elements 

Again  neglecting  source  terms,  the  equations  of  motion  for  elastic  medium  for  the  displacement 
components  u  and  v  are 

u'^pu  +  (ciiu,!  +  ci2U,z)^^  +  C33  =  0  (35a) 

J^pv  +  C33  (u  z  +  v^x)^x  +  +  C22V  2)  ^  =  0  (35b) 


where  the  density  p  and  the  elastic  constants  cn  =  C22  =  A+2/x,  C12  =  A  and  C33  =  p  are  all  spatially 
varying.  The  subscripts  and  denote  partial  derivatives  with  respect  to  x  or  respectively. 

As  in  the  acoustic  case  (8),  the  displacements  inside  the  elements  ^  are  approximated  by 
interpolation  functions  (Zienkiewicz,  1977;  Schwarz,  1988;  Murphy  and  Chin-Bing,  1991): 


4 


■“  =  E  W 

j  =  l 

4 

(36a) 

V  =  Y^VjNj{x) 

(36b) 

The  complex  valued  weighting  coefficients  Uj  and  Uj  are  the  components  of  the  displacements  at 
the  element’s  corners  x^.  The  shape  functions  iVj(x)  are  the  same  as  in  the  acoustic  case,  e.g. 
iV3(x)  =  xz  (Figure  2). 

In  the  equation  of  motion  (35),  we  replace  the  displacements  u  and  v  by  the  interpolations  (36). 
Applying  Galerkin’s  method  for  each  element  fl,  we  multiply  the  resulting  expression  by  the  test 
function  Ni{x)  and  integrate  over  the  element  Q. 


'h  dA  JJ ^  NiNj  Uj  + 

E  {//j^  +  C33^^^^ 


xcrn  dl  =  0  (37a) 


,x^j,x  +  C22Ni^zNj^z^  dA  —  j j .  NiNj  vj  ~  J,  Ni 


zcrn  dl  =  0  (37b) 
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where  we  used  the  divergence  theorem  to  transform  some  of  the  volume  integrals  into  line  integrals. 
The  quantity  cr  denotes  the  stress  tensor  along  the  boundary.  Thus,  evaluating  (37)  for  all  iVj(x) 
will  yield  a  set  of  equations  for  the  unknown  node  variables  Uj  and  vj.  Again,  equation  (37)  defines 
the  stiffness  matrices  S*-’,  mass  matrices  M”  and  the  force  vectors  T: 


=  Mlf=  Jj ^pu^  NiNjdA 

(38) 

—  JJ ~  X  C33ffi,zf^j,z'j  dA 

(39) 

q22 

^ij 

=  jj ^  {czsNi^xNj^x  +  C22Ni,zNj,zj  dA 

(40) 

q12 

^ij 

(41) 

fl 

—  j  Nixcrn  dl 

(42) 

fi 

=  J_  Nizcrn  dl 

(43) 

If  the  density  p  and  the  elastic  constants  Cn,  C22,  C12  and  C33  are  treated  as  constants  within 
each  element  fi,  all  integral  (38-41)  can  be  evaluated  exactly.  For  elements  which  are  in  the 
interior,  the  boundary  integrals  (42,43)  are  not  zero,  but  its  contributions  will  exactly  cancel  with 
like  terms  coming  from  neighboring  elements  because  displacements  and  stresses  are  continuous 
across  boundaries.  Therefore,  the  line  integrals  have  to  be  taken  into  account  only  on  the  domain 
boundary.  Assuming  that  xah  and  zcrh  can  be  approximated  by  functions  similar  to  Ni  along  the 
boundary,  we  can  replace  (42,42)  by 

fi  =  ^FijX(r{-Kj)n  _  (44) 

i=i 

fi  =  Y^FijZO-{yij)h  (45) 

jzzl 

Fij  =  J^NiNjdl  (46) 

which  can' also  be  evaluated  analytically.  Mapping  all  the  local  contributions  of  M“  and  f* 
into  global  node  numbers  yields  the  global  matrices  M”,  f*  and  the  global  nodal  vectors  u  and 
V.  Writing  the  global  matrix  system  reduces  to  the  simpler  system 

-I- _  fi  ^  0 

K^^u-t-K^V-f^  =0. 
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(47a) 

(47b) 


Coupling  the  Regions 

Each  vector  u  and  v  obtained  from  the  finite  elements  is  split  into  two  subvectors  u^,  and 
v^,  v^,  respectively.  The  node  variables  Ui  and  vt  from  inside  the  domain  are  collected  in  the 
vectors  and  v^,  respectively.  The  subvectors  and  accommodate  the  node  variables  whose 
nodal  points  Xi  lie  on  the  boundary  Since  the  boundary  belongs  to  both  domains,  the 

respective  wavefields  w®  and  have  to  match  along  the  boundary  and  we  can  replace  the  node 
variables  uf  and  vf  by 

j° 

uf  =  -Y,  {hut (xj )  +  i>kUk  (xj ) }  +  (xj )  (48a) 

k=l 

J° 

vf  =  z-Y  [hVk  (Xj)  +  ipkVk  (Xj)}  +  v’"‘'(xj).  (48b) 

k=l 

Furthermore,  we  find  cr(xj)  by  evaluating 


=  Y  {h<rt (Xj)  +  i’ktrf  (Xj)}  +  o-’”‘=(xj)  (49) 

fc=l 

Combining  (47),  (48)  and  (49)  yields  the  coupled  MMP-FEM  system.  For  the  sake  of  clarity,  we 
explicitly  expand  (47a): 


where  is  the  number  of  node  variables  inside  the  heterogeneous  region  is  the  node 

number  of  the  first  nodal  point  lying  on  the  boundary  d^l^.  The  value  of  is  + 1.  As  before,  J 
is  the  total  number  of  node  points.  Finally,  J®  is  the  total  number  of  functions  used  for  the  MMP 
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expansion  of  the  outside  field.  The  rows  of  (47b)  follow  the  same  outline.  The  resulting  combined 
system  of  equations  is  of  similar  form  as  (21)  and  can  be  solved  using  the  same  scheme. 


(  K{i  K{2  \  /  u  \  /  ff 

KI2  ^  ^  w  ^  ii 

Kfi  Kf2  -J-f  ‘  <^  fP 

\kO  Kg  \^fO 

As  in  the  acoustic  case,  the  submatrices  resulting  from  the  interior  problem  are  sparse  and 
square.  All  other  submatrices  are  rectangular.  We  reduce  the  above  system  by  Gaussian  eliinination 
of  the  node  variables  u  and  v  which  yields  a  new  system. 


(52) 


The  lower  half  of  (52)  can  now  be  solved  in  the  least  squares  sense  by  QR  decomposition. 


(53) 


The  node  variables  in  the  heterogeneous,  interior  region  are  recovered  by  back-substitution,  the 
upper  half  of  (52)  is  already  in  upper  triangular  form. 


K{2  \  /  u  \  \ 

0  Kg  )  \v  )  / 


IMPLEMENTATION 


Because  the  technique  is  a  mixture  of  MMP  expansions  and  the  FE  method,  the  finite  element 
method  is  merged  into  the  prior  MMP  codes  (Imhof,  1995a, b).  Thus,  the  method  is  implemented 
on  a  nCUBE2  parallel  computer  using  the  computer  language  C++.  The  object  oriented  design  has 
the  advantage,  that  the  coupling  as  described  in  (20)  and  (50)  is  batsically  hidden  in  objects  for  node 
variables,  finite  elements  and  the  expansion  functions  for  the  exterior.  First,  the  objects  for  the  finite 
elements  calculate  the  local  M,  S  and  F  matrices.  Then,  the  resulting  coefficients  have  mapped 
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into  the  global  equation  system.  Objects  for  internal  node  variables  simply  map  the  coefficients 

for  the  pj  into  the  global  system  of  equations.  In  contrast,  objects  for  node  variables  on  the 
boundary  automatically  evaluate  the  MMP  expansion  at  the  node  point  as  described  in  equations 
(18), (19)  or  (48), (49),  weight  the  expansion  with  the  appropriate  or  coefficient  and  map 
the  resulting  coefficients  for  or  into  the  global  system. 

To  reduce  numerical  noise,  the  materials  are  made  slightly  lossy  by  adding  a  small  imaginary 
component  to  the  angular  frequency.  If  seismograms  are  calculated  by  Fourier  synthesis,  the 
true  amplitude  is  recovered  by  a  multiplication  with 

NUMERICAL  RESULTS:  ACOUSTICS 

As  a  first  test,  we  simply  embed  a  homogeneous  region  in  a  homogeneous  fullspace.  The  wavefield 
in  the  embedded  region  is  modelled  by  FE.  The  wavefields  in  the  fullspace  are  expanded  into  a 
MMP  series.  The  material  parameters  in  both  regions  are  the  same.  Hence,  all  coefficients  of  the 
MMP  expansion  should  be  zero,  while  the  FE  solution  should  simply  interpolate  the  incoming  field. 
Clearly,  due  to  the  discretization  of  the  field  in  the  interior,  the  solution  in  the  interior  will  deviate 
from  the  incident  field  and  thus,  an  additional  scattered  field  will  be  induced.  The  strength  of 
this  induced  field  is  both  a  function  of  the  number  of  elements  per  wavelength  and  the  angle  of 
incidence  of  the  source  field.  The  embedded  region  consists  of  18  *  18  elements,  each  4m  *  4m  in 
size.  The  MMP  expansion  is  (2)  with  M  =  4  and  N  =  4.  Altogether,  36  expansion  functions  are 
used.  Figure  3  shows  the  exact  position  of  node  points  and  expansion  centers.  The  source  field  is 
a  plane  wave,  where  the  angle  of  incidence  ranges  from  0°  up  to  45°.  As  a  measure  for  the  error, 
we  use  <  (P  -  p*"<:)/P'"<=)  >  al&ng  the-bpundary  of  the  inclusion.  Starting  with  250  elements  per 
wavelength  (EPW),  the  number  of  EPW  is  steadily  decreased  down  to  2  EPW.  The  results  are 
shown  in  Figure  4.  As  can  be  seen,  the  error  increases  slowly  until  the  induced  fields  are  of  similar 
size  to  the  source  field.  Using  10  EPW  yields  an  error  of  about  5%.  Also,  the  error  becomes  slightly 
smaller  the  more  the  incident  plane  wave  propagates  in  the  diagonal  direction. 

To  test  the  accuracy  of  the  MMP-FE  technique,  we  compare  the  scattering  from  an  acoustic 
cylinder  with  the  well-known  analytical  series  solution  (Pao  and  Mow,  1973).  The  velocity  inside 
the  cylinder  is  3000m/s;  the  velocity  outside  the  cylinder  the  velocity  is  2000m/s.  In  both  regions, 
the  density  is  kept  constant  at  2000kg/m^.  The  radius  a  of  the  cylinder  is  44m.  To  simplify  the 
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generation  of  the  FE  rhesh,  a  square  region  larger  than  the  actual  cylinder  is  discretized  by  24 
elements  in  either  direction.  Due  to  the  symmetry  of  the  problem,  only  one  multipole  (2),  where 
M  =  1  and  N  =  20,  located  at  the  origin,  is  used.  The  geometry  is  shown  in  Figure  5.  The  size 
of  the  elements  is  4m.  Two  different  wavelengths  are  used:  25m  and  100m.  Magnitude  and  phase 
are  presented  in  Figures  6  and  7.  In  the  case  of  ka  =  10,  the  deviations  of  the  FE-MMP  solution 
from  the  analytical  one  are  due  to  the  finite  element  size.  Reducing  the  element  size  reduces  the 
deviations.  Furthermore,  the  largest  deviations  correlate  with  the  smallest  magnitudes  as  can  be 
observed  in  both  Figures  6  and  7.  This  is  an  effect  of  the  least-squares  solving  procedure.  The 
solver  uniformly  minimizes  the  misfit  at  each  boundary  point.  Thus,  if  the  average  misfit  is  e,  any 
true  field  value  smaller  than  e  is  lost  in  the  misfit.  If  better  accuracy  is  desired,  the  solution  should 
be  calculated  again  with  the  equations  scaled  by  the  reciprocal  field  obtained  before.  Basically, 
and  should  be  scaled  by  Further  details  on  scaling  can  be  found  in  the  prior  paper 
(Imhof,  1995a). 

Lastly,  we  calculate  the  seismogram  for  a  complex  geometry  depicted  in  Figure  8.  The  scatterers 
are  roughly  180m  long  and  35m  thick.  The  velocity  and  density  in  the  background  are  2000m/s 
and  2000kg/m^,  respectively.  The  velocity  and  the  density  in  the  two  scatterers  are  3000m/s  and 
2000kg/m^,  respectively.  Each  finite  element  is  3m  by  3m  in  size.  For  each  scatterer,  five  centers 
of  expansion  are  used.  At  each  center  Xm,  expansion  of  the  form 

n=:*-8 

is  set  up.  The  incident  field  p  is  an  explosive  line  source  modulated  with  a  Ricker  wavelet 
(Ricker,  1977)  of  50  Hz  center  frequency.  Altogether,  64  receivers  will  measure  the  pressure  of  the 
scattered  field.  The  resulting  seismogram  is  shown  in  Figure  9. 

NUMERICAL  RESULTS:  ELASTICS 

As  in  the  acoustic  case,  the  first  test  is  to  embed  a  homogeneous  region  in  a  homogeneous  fuUspace. 
The  wavefields  in  the  embedded  region  are  modelled  by  FE,  while  the  wavefields  in  the  fullspace 
are  expanded  into  a  MMP  series.  Because  the  material  parameters  in  both  regions  are  the  same,  all 
coefficients  of  the  MMP  expansion  should  be  zero  and  the  FE  solution  should  perfectly  interpolate 
the  incoming  field.  Clearly,  due  to  the  discretization  of  the  fields  in  the  interior,  the  solution  will 
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deviate  from  the  incident  field  and  thus,  additional  scattered  fields  will  be  induced.  The  strength 
of  these  induced  fields  is  both  a  function  of  the  number  of  elements  per  wavelength  and  the  angle 
of  incidence  of  the  source  field.  The  embedded  region  consists  of  18  *  18  square  elements,  each 
4m  *  4m  in  size.  The  MMP  expansion  is  the  same  as  (29)  with  M  =  4  and  =  4.  Altogether, 
2*36  expansion  functions  are  used.  Figure  3  shows  the  exact  position  of  node  points  and  expansion 
centers. 

As  source  fields,  we  use  plane  waves  of  purely  P  or  S  polarization.  Either  experiment  is  per¬ 
formed  twice,  first  for  an  angle  of  incidence  of  0°,  then  45°.  Scanning  through  a  range  of  frequencies 
allows  is  to  see  how  polarization,  orientation  of  the  elements,  and  the  number  of  elements  per  wave¬ 
length  affect  the  solution.  As  a  measure  for  the  error,  we  use  <  (|w  —  w*”°|)/|w*”°|  >  along  the 
boundary  of  the  inclusion.  Again,  we  start  with  250  elements  per  P-wavelength  (EPW)  and  de¬ 
crease  the  number  of  EPW  down  to  2.  The  results  are  shown  in  Figure  10.  As  expected,  the  errors 
increase  with  increasing  frequency.  Interestingly,  incident  S  waves  of  very  low  frequency  are  less 
affected  than  incident  P  waves  of  the  same  frequency.  But  with  increasing  frequency,  the  rate  with 
which  the  error  grows  is  larger  for  incident  S-  than  incident  P-waves.  If  less  than  5  EPW  are  used, 
the  S-waves  are  aliased  and  the  results  become  meaningless.  In  general,  waves  incident  at  45°  are 
less  affected  by  the  grid  size  than  waves  incident  in  the  normal  direction.  Using  25  EPW  yields  an 
error  of  about  8%. 

To  test  the  accuracy  of  the  MMP-FE  technique  in  the  elastic  case,  we  compare  the  scattering 
from  a  cylinder  with  the  analytical  series  solution  (Pao  and  Mow,  1973).  Inside  the  cylinder,  the 
P-velocity  is  3000m/s  and  the  S-velocity  is  1700m/s.  Outside  the  cylinder  the  P-velocity  is  2000m/s 
and  the  S-velocity  is  I300m/s.  In  both  regions,  the  density  is  2000kg/m^  and  the  Poisson’s  ratio 
is  The  radius  a  of  the  cylinder  is  12m.  To  simplify  the  generation  of  the  FE  mesh,  a  square 
region  larger  than  the  actual  cylinder  is  discretized  by  24  elements  in  either  direction.  Due  to  the 
symmetry  of  the  problem,  only  one  multipole  (29),  where  M  =  \  and  N  =  20,  located  at  the  origin, 
is  used.  The  geometry  is  shown  in  Figure  5.  The  size  of  the  elements  is  Im.  The  wavelength  of 
the  incident  P-wave  is  50m  and  the  wavelength  of  the  incident  S-wave  is  32m.  The  magnitude  and 
phase  of  the  u  and  the  v  components  are  shown  in  Figures  11  and  12.  For  all  incident  phases, 
the  match  between  the  analytical  solution  and  the  results  obtained  from  the  MMP-FE  method  are 
excellent. 
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Lastly,  we  calculate  the  seismogram  for  a  complex  geometry  depicted  in  Figure  13.  A  scatterer 
is  illuminated  by  a  line  source.  The  scatterer  is  roughly  180m  long  and  35m  thick.  The  P-  and 
S-velocities  and  density  in  the  background  are  2000m/s,  1300m/s  and  2000kg/m^,  respectively.  The 
P-  and  S-velocities  and  density  in  the  scatterer  are  3000m/s,  1730m/s  and  2000kg/m^,  respectively. 
The  Poisson’s  ratio  is  cr  =  0.25  in  both  the  scatterer  and  in  the  background.  Each  finite  element 
is  2m  by  2m  in  size.  In  the  scatterer,  five  expansions  of  the  form  (29)  with  M  —  5  and  N  =  7  are 
used.  Two  different  incident  fields  are  chosen:  a  compressional  and  a  rotational  line  source.  Each 
source  is  modulated  with  a  Ricker  wavelet  (Ricker,  1977)  of  50  Hz  center  frequency.  Altogether,  64 
receivers  measure  the  vertical  displacement  component  of  the  total  field.  The  resulting  seismograms 
are  shown  in  Figures  14  and  15. 


SUMMARY 

The  MMP  code  has  been  successfully  coupled  with  the  FE  method  in  both  acoustic  and  elastic 
media.  The  coupling  of  the  two  methods  enhances  their  usefulness  for  a  range  of  problems.  The 
FE  technique  allows  the  simulation  of  wave  propagation  in  heterogeneous  materials.  The  MMP 
expansions  allow  to  calculate  propagating  waves  in  homogeneous  (unbounded)  regions  in  an  efiicient 
manner  because  they  commonly  need  less  unknowns  to  be  evaluated  and  solved  for  than  comparable 
methods. 

Steady-state  solutions,  as  well  as  seismograms  obtained  by  Fourier  synthesis,  were  calculated 
for  a  range  of  different  problems  for  both  acoustic  and  elastic  media.  Where  available,  the  solutions 
obtained  by  the  combined  MMP-FEM  scheme  compared  favorably  with  the  analytical  solutions. 

The  combined  scheme  compensates  for  the  individual  weaknesses  of  MMP  and  FEM  and  takes 
advantage  of  both  their  strengths.  Thus,  the  method  is  well-suited  to  solve  two-dimensional  scat¬ 
tering  problems  for  a  range  of  problems  which  neither  method  could  handle  alone. 
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Figure  1;  The  generic  scattering  problem  to  be  solved  by  the  hybrid  MMP-FEM  technique.  A 
heterogeneous  scatterer  is  embedded  in  a,  homogeneous  background  Q,^.  In  the  acoustic  case, 
the  incident  field  is  and  the  scattered  field  is  p^.  In  the  elastic  case,  the  incident  field  is 
and  the  scattered  field  The  triangles  symbolize  expansion  centers  for  the  MMP. 
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Figure  2:  The  shape  function  Nz{x)  =  xz  for  a  square  unit  element  Cl  and  bilinear  interpolation. 
The  other  shape  functions  iVi(x),  N2{x.)  and  iV4(x)  are  obtained  by  rotations  of  -180°,  -90°  and 
90°,  respectively. 
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Figure  3:  An  embedded  homogeneous  region.  Each  dot  represents  a  node  point  and  each  di¬ 
amond  a  MMP  expansion  center.  From  each  expansion  center  we  set  up  an  expansion 
EL-4  (ko  I*  -  x^l). 
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Boundary  Error  as  a  Function  of  Element  Size  and  Orientation 


Figure  4:  Relative  boundary  error  <  (p  — >  as  a  function  of  the  number  of  elements  per 
wavelength  (EPW)  and  propagation  angle  of  the  incident  field  with  respect  to  the  finite  elements. 
10  EPW  yields  an  error  of  5%.  Waves  propagating  diagonally  are  less  affected  by  larger  element 
sizes. 
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Figure  5:  A  cylindrical  scatterer  is  illuminated  by  a  plane  wave.  Outside,  the  velocity  is  2000m/s; 
inside,  the  velocity  is  3000m/s.  The  grid  represents  the  finite  elements  used.  The  grid  spacing  is  4m 
and  the  radius  of  the  cylinder  is  44m.  The  triangle  denotes  the  expansion  center  for  the  multipole. 


Cylindrical  Scatterer:  Magnitude 


Figure  6:  A  cylindrical  scatterer  is  illuminated  by  a  plane  wave  propagating  in  the  positive,  hori¬ 
zontal  direction.  Shown  is  the  magnitude  |P|  for  ka  =  2.5  and  fca  =  10  as  a  function  of  angle. 
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Cylindrical  Scatterer:  Phase 


Figure  7:  A  cylindrical  scatterer  is  illuminated  by  a  plane  wave  propagating  in  the  positive,  hori¬ 
zontal  direction.  Shown  is  the  phase  arg(P)  for  ka  =  2.5  and  ka  =  10  as  a  function  of  angle. 
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Figure  8:  Generic  scatterers  used  to  calculate  seismograms.  Two  homogeneous  scatterers  are 
embedded  an  a  homogeneous  background.  The  velocity  and  the  density  in  the  background  are 
2000m/s  and  2000kg/m^,  respectively.  The  velocity  and  the  density  in  the  two  scatterers  are 
3000m/s  and  2000kg/m^,  respectively.  The  triangles  show  the  location  of  the  centers  for  the  MMP 
expansion. 
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Boundary  Error  as  a  Function  of  Element  Size,  Orientation  and  Type  of  Incident  Field 


Figure  lOr  Relative  boundary  error  <  |w  —  >  as  a  function  of  the  number  of  elements 

per  P-wavelength  (EPW)  and  propagation  angle  of  the  type  of  incident  field  with  respect  to  the 
finite  elements.  25  EPW  yields  an  error  of  8%.  Waves  propagating  diagonally  are  less  affected 
by  larger  element  sizes.  Also,  an  incident  S  wave  is  more  affected  by  the  element  size  because  its 
wavelength  is  roughly  half  as  long  as  the  corresponding  P  wave’s. 
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Figure  11:  A  cylindrical  scatterer  is  illuminated  by  a  plane  wave  propagating  in  the  positive, 
horizontal  direction.  Both  P  and  S  modes  are  used  as  incident  fields.  Shown  are  the  magnitude  of 
\u\  and  |t;|  for  ka  =  1.5  and  la  =  2.3  as  a  function  of  angle. 


Cylindrical  Scatterer:  Phase 


Figure  12:  A  cylindrical  scatterer  is  illuminated  by  a  plane  wave  propagating  in  the  positive, 
horizontal  direction.  Both  P  and  S  modes  are  used  as  incident  fields.  Shown  are  the  magnitude  of 
arg(w)  and  arg(t;)  for  fca  =  1.5  and  la  =  2.3  as  a  function  of  angle. 


87 


I  I  I  I  I  I  I  1 1 1 1  I  I  1 1 1  I  I  M  I  I  I 

A 


H  I  I  M  I  l)Kl  M  I  1 1 1  I  1 1  I  I  I  1 1 1  I  I  I  I  I  I  I  1 1  I  I  I  I  I  I 


compressional  /  rotational 
line  source 

(x ,  co) 


64  receivers  ->  <r 
5m 


100m 


- 

>C'' 

St 

ii 

% 

ii 

■  V - 

m 

X 

•  j 

s 

Si 

1” 

3 

— 

ij 

< - 

180m 


A 

35m 

Y 


Figure  13;  Generic  scatterer  used  to  calculate  elastic  seismograms.  The  scatterer  is  embedded  in 
a  homogeneous  background.  The  P-  and  S- velocities  and  density  in  the  background  are  2000m/s, 
1300m/s  and  2000kg/m^,  respectively.  The  P-  and  S-velocities  and  density  in  the  scatterer  are 
3000m/s,  1730m/s  and  2000kg/m^,  respectively.  Thus,  Poison’s  ratio  in  the  scatter  is  the  same  as 
in  the  background.  The  triangles  show  the  location  of  the  centers  for  the  MMP  expansion. 
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Figure  14:  The  vertical  component  of  the  seismogram  for  the  complex  geometry  depicted  in  Fig¬ 
ure  13.  The  incident  field  is  a  compressional  line  source. 
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Figure  15:  The  vertical  component  of  the  seismogram  for  the  complex  geometry  depicted  i 
13.  The  incident  field  is  a  rotational  line  source. 
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ABSTRACT 

In  this  paper  we  compared  scattered  waves  from  2-D  and  3-D  interface  structures.  The  modeling 
technique  is  the  3-D  time  domain  finite  difference  method.  The  scheme  is  second-order  accurate 
in  time  and  fourth-order  accurate  in  space.  It  is  implemented  on  a  massively  parallel  nCUBE 
computer.  In  order  to  investigate  the  characteristics  of  2-D  and  3-D  rough  surface  scattering,  we 
consider  an  acoustic-elastic  boundary,  which  is  described  by  a  Gaussian  autocorrelation  function. 
The  F-K  analysis  of  reflected  signals  shows  that  2-D  scattering  generates  similar  amounts  of  forward 
and  back  scattering,  while  in  the  3-D  case,  more  forward  and  less  back  scattering.  The  3-D  effects 
also  show  larger  reflected  energy  than  the  2-D  case,  especially  near  the  normal  incident.  The 
out-of-plane  scatterings  are  clearly  demonstrated  on  the  F-K  spectra  in  the  3-D  case.  In  the  2-D 
simulations,  we  have  to  keep  in  mind  that  it  tends  to  overestimate  the  amount  of  backscattering 
energy. 
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INTRODUCTION 


Many  boundaries  in  the  Earth  may  be  characterized  as  irregular  interfaces,  such  as  the  ocean  floor, 
the  Moho  discontinuity,  and  surface  terrain.  The  wave  scattering  from  interfaces  has  become  a 
subject  of  wide  interest.  The  effects  of  scattering  on  the  observations  can  often  be  identified  on 
amplitudes,  traveltimes,  spectra,  and  waveforms.  Although  scattering  from  smooth  interfaces  is 
easily  modeled  and  well-understood,  the  roughness  inherent  in  many  natural  interfaces  introduces 
an  incoherent  component  into  the  wavefield,  which  is  poorly  understood.  A  number  of  theoretical 
approaches  have  been  developed  to  help  us  understand  the  rough  interface  scattering  problem. 
Reviews  of  these  works  were  given  by  Ogilvy  (1987)  and  Thorsos  and  Jackson  (1989).  These 
studies  are  based  on  the  perturbation  method.  They  require  interface  properties,  such  as  the 
overall  topography  (amplitude  and  slope),  to  be  small  compared  with  the  incident  wavelength. 
Under  the  approximation  these  studies  give  analytical  solutions  to  the  problems.  Full-waveform 
techniques,  such  as  the  boundary  integral  method,  are  implemented  to  model  completely  the  wave 
interaction  with  the  interface  topography  (Bouchon,  1985;  Kawase,  1988;  Gerstoft  and  Schmidt, 
1991;  Sanchez-Sesma  and  Campillo,  1991).  The  finite  difference  method  has  also  become  effective 
in  studying  scattering  from  a  rough  interface  (Lavender  and  Hill,  1985;  Burns  and  Stephen,  1990; 
Dougherty  and  Stephen,  1991;  Fricke,  1993). 

Most  numerical  modelings  have  focused  on  two-dimensional  models,  such  as  the  ocean  floor,  ice 
sheet,  and  shallow  crust.  These  2-D  simulations  have  demonstrated  the  basic  effects  of  scattering 
and  wave  conversions  at  interfaces.  Nature  boundaries,  however,  are  three  dimensional.  The 
3-D  scattering  features,  such  as  the  out-of-plane  effect,  cannot  be  addressed  adequately  by  2-D 
modeling.  It  is  important  to  understand  the  difference  between  2-D  and  3-D  surface  scattering. 
This  knowledge  can  help  us  assess  the  limitations  of  2-D  modeling.  It  is  the  purpose  of  this  paper 
to  show  an  example  of  the  differences  between  2-D  and  3-D  rough  surface  scattering.  We  chose 
the  finite  difference  method  as  the  modeling  technique.  The  finite  difference  method  can  propagate 
complete  wave  fields  through  arbitrarily  complex  media.  However,  the  finite  difference  method  is 
very  computationally  intensive,  especially  for  the  3-D  problem.  Parallel  computing  makes  the  finite 
difference  method  a  very  attractive  choice.  The  massively  parallel  computer  is  used  in  this  paper 
to  implement  the  3-D  finite  difference  scheme.  We  chose  the  rough  interface  as  an  acoustic-elastic 
boundary.  The  interface  is  changed  from  1-D  (flat)  to  2-D  (variation  in  one  direction)  and  to  3-D 
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(variations  in  two  directions)  to  do  the  numerical  simulations.  This  enables  us  to  compare  the 
scattering  from  2-D  and  3-D  rough  surfaces  directly. 


3-D  FINITE  DIFFERENCE  METHOD 


where  v  is  velocity,  p  is  density,  and  r  is  stress.  A  and  p  are  Lame  constants.  The  reason  for 
formulating  the  second-order  wave  equations  into  the  first-order  hyperbolic  system  of  equations  is 
that  once  this  system  is  discretized  on  a  staggered  grid,  it  is  valid  for  any  Poisson’s  ratio  (Virieux, 
1986).  The  fluid-solid  boundary  can  be  treated  simply  by  setting  shear  modulus  to  zero. 

The  first-order  hyperbolic  equations  (1)  and  (2)  are  discretized  on  a  staggered  grid  (Madariaga, 
1976),  which  is  shown  in  Figure  l.We  approximate  the  first-order  time  derivative  by  second  order 
centered  finite  difference  and  the  first-order  spatial  derivatives  by  fourth-order  centered  finite  dif¬ 
ference.  The  media  parameters  p,  X  and  p  are  given  at  the  grid  point  where  the  normal  stresses 
Txx,Tyy,Tzz  are  assigned  (see  Figure  l).In  the  calculation  to  update  velocities,  the  needed  density 
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values  are  obtained  from  the  average  of  the  two  assigned  densities  nearby.  The  shear  moduli  used 
to  update  the  shear  stress  are  determined  by  the  harmonic  average  of  the  four  shear  moduli  nearby 
instead  of  the  arithmetic  average  (Kostek,  1991).  This  harmonic  average  method  can  automatically 
put  the  shear  modulus  zero  at  the  fluid-solid  boundary. 

The  dispersion  analysis  shows  that  the  fourth-order  flnite  diflFerence  has  much  less  dispersion 
and  grid  anisotropy  than  the  second-order  one  fot  both  P  and  S  waves.  The  rule  of  thumb  is  that 
we  need  five  samples  per  wavelength  for  the  fourth-order  finite  difference  to  control  the  dispersion 
and  the  anisotropy  at  less  than  1%.  The  stability  condition  of  the  scheme  is  given  by  (Cheng,  1994) 


under  the  condition  of  grid  size  Ax  =  Ay  =  Az  =  A,  where  tji  =  g  and  ^2  =  ^'^e  the 

coefficients  of  the  fourth-order  finite  difference  approximation  to  the  first-order  derivative,  a  is 
the  maximum  P  wave  velocity  in  the  model.  The  absorbing  boundary  condition  is  applied  to  the 
outside  boundaries  of  the  grid  to  minimize  the  reflections.  Higdon’s  absorbing  boundary  condition 
is  used  (Higdon,  1990). 

Parallel  computing  provides  a  new  means  to  overcome  the  limitations  of  the  memory  and  the 
speed  of  a  single  processor.  In  the  finite  difference  method  all  the  calculations  involve  only  local 
interactions  of  the  velocities  and  stresses.  For  example,  in  the  fourth-order  finite  difference  scheme 
only  two  nearby  grid  points  data  are  needed  to  update  the  current  grid  point.  This  can  be  eflficiently 
executed  on  a  multiple  instruction  and  multiple  data  (MIMD)  parallel  computer.  A  staggered-grid 
fourth-order  finite  difference  scheme  discussed  above  is  paralleled  on  nCUBE  2S  massively  parallel 
computer  (Cheng,  1994). 


NUMERICAL  EXAMPLES  AND  F-K  ANALYSIS 


Test  of  the  Finite  Difference  Method 

We  consider  a  fluid-filled  cylindrical  borehole  embedded  in  a  homogeneous  elastic  formation.  The 
finite  difference  results  are  compared  with  the  discrete  wavenumber  method.  The  physical  property 
of  the  fluid  and  solid  are  listed  in  Table  1.  The  borehole  radius  is  0.1  m.  A  point  explosion  source 
and  the  pressure  receivers  are  located  at  the  center  of  the  borehole.  The  Kelly  wavelet  at  center 
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frequency  14  kHz  is  used  as  the  source  time  function  (Kelly  et  al^  1976).  A  grid  size  of  0.01  m  and 
a  time  step  of  0.001  ms  were  chosen.  The  second-order  Higdon’s  absorbing  boundary  condition  is 
used  to  reduce  the  artificial  boundary  reflections.  The  local  P  and  S  wave  velocities  are  used  at 
the  absorbing  boundary.  A  grid  of  70  x  70  x  200  is  used  to  build  the  borehole  model.  The  finite 
difference  results  are  compared  with  the  discrete  wavenumber  solutions  (Cheng  and  Toksoz,  1981*) 
in  Figure  2.There  are  body  waves  (refracted  P  and  S)  and  guided  waves  (Stoneley  and  pseudo- 
Rayleigh)  propagated  inside  the  borehole.  The  comparison  shows  excellent  agreement  between  the 
two  solutions.  This  test  shows  the  finite  difference  scheme  handled  the  sharp  fluid-solid  boundary 
very  well.  Even  if  the  borehole  is  approximated  by  a  rough  edged  circle  in  the  finite  difference 
method,  we  still  can  get  a  good  match  between  the  two  methods. 

Synthetic  Reflection  Data 

In  order  to  investigate  the  reflection  characteristics  from  the  rough  surface,  we  consider  an  irregular 
fluid-solid  interface.  There  are  only  reflected  P  waves  in  the  fluid.  This  will  greatly  simplify  the 
analysis  of  reflected  signals.  The  property  of  the  fluid  and  the  solid  are  listed  in  Table  1  and  a  3-D 
view  of  the  model  is  shown  in  Figure  3. The  source  time  function  of  a  point  explosion  is  a  Kelly 
wavelet  with  a  center  frequency  of  15  Hz.  We  choose  the  grid  discretization  interval  to  be  one-tenth 
the  shortest  wavelength  at  the  source  center  frequency,  which  is  100  m.  The  time  step  is  0.001  s 
and  the  size  of  the  model  is  200  x  200  x  100. 

In  this  study,  we  place  a  2-D  pressure  receiver  array  (64  x  64)  in  a  plane  parallel  to  the  interface 
and  on  one  side  of  the  source.  The  array  receiver  spacing  is  20  m.  We  adopt  a  Gaussian  surface  to 
describe  the  fluid  and  solid  boundary.  It  has  a  Gaussian  spatial  correlation  function  and  a  Gaussian 
distribution  about  the  mean  (Frankel  and  Clayton,  1986).  This  interface  can  be  described  by  two 
properties:  the  correlation  length,  a,  and  the  standard  deviation  in  height,  6.  The  correlation 
length  corresponds  to  the  average  distance  between  nearby  peaks  and  valleys  along  the  interface, 
and  the  standard  deviation  gives  the  root  mean  square  (rms)  deviation  of  the  interface  height  from 
its  mean,  which  is  set  to  zero.  The  rms  slope  <f)  of  the  interface  is  given  by 

c 

tan(<;!!i)  =  v^-  (4) 

0/ 

where  the  correlation  length  is  set  equal  to  the  central  incident  wavelength  (100  m)  and  rms  slope 
is  set  at  30 
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A  3-D  rough  surface  is  introduced  with  the  correlation  length  in  the  X  and  Y  direction  both 
equal  to  100  m,  while  the  2-D  rough  interface  has  a  correlation  length  in  the  X  direction  of  100  m, 
and  of  oo  in  the  Y  direction.  The  2-D  Gaussian  rough  surface  is  plotted  in  Figure  4b. and  the  3-D 
rough  surface  is  plotted  in  Figure  4b.The  2-D  rough  surface  is  inherited  from  the  3-D  one.  It  is 
simply  obtained  by  extending  the  3-D  rough  surface  at  Y=1  in  the  Y  direction.  This  will  make  the 
amplitude  comparison  interesting  later.  In  that  case  we  obtain  the  amplitudes  from  3-D  scattering 
at  the  plane  which  has  the  exactly  same  surface  profile  as  the  2-D  one.  The  source-receiver  plane 
is  600  m  above  the  mean  height  of  the  rough  interface. 

We  have  calculated  the  synthetic  data  set  first  for  a  flat  interface,  as  shown  in  Figure  5a.Figures 
5b  and  c  show  the  seismogram  obtained  from  2-D  and  3-D  rough  interface  models.  These  plots  are 
a  2-D  slice  of  the  3-D  data  set  for  fixed  Y.  The  direct  arrivals  from  the  finite  difference  calculation 
are  removed,  and  they  are  not  shown  in  the  plots.  The  coda  waves  in  Figures  5b  and  c  show  the 
effect  of  the  rough  interface.  The  coda  generated  from  the  2-D  rough  interface  have  more  structures 
than  the  3-D  one.  Also  the  clear  reflections  from  the  flat  interface  are  smeared  in  the  2-D  and  3-D 
rough  surface  cases.  The  F-K  analysis  of  the  reflected  signals  will  follow. 

F-K  Analysis  of  the  Reflected  Data 

F-K  analysis  is  a  useful  technique  to  illustrate  the  magnitude  and  direction  of  energy  arriving  at 
the  array.  It  is  widely  used  in  the  seismic  coda  analysis,  and  Dainty  and  Toksoz  (1990)  give  a  good 
summary  of  this  approach.  In  this  case  we  use  a  simple  3-D  Fourier  transform  algorithm  on  the 
synthetic  data,  transforming  in  both  spatial  directions,  and  in  time.  The  data  is  displayed  at  a 
given  frequency  to  show  the  magnitude  and  direction  of  the  scattered  waves.  The  source  center 
frequency  will  be  of  interest  here.  Before  we  analyze  the  reflected  data  from  the  interfaces,  we  first 
look  at  the  direct  arrivals  (Figure  6a). The  direct  arrivals  are  propagated  straight  from  the  source 
to  the  receivers  on  the  same  plane.  A  radius  of  10  corresponds  to  the  direct  wave  arriving  with  a 
horizontal  phase  velocity  of  1.5  km/s  at  15  Hz.  Only  half  of  the  circle  is  shown  on  the  figure.  This 
is  because  the  receivers  are  only  on  one  side  of  the  source.  For  the  purpose  of  comparison  the  F-K 
spectra  of  the  reflections  from  the  flat  surface  is  shown  in  Figure  6b.  The  small  wavenumber  range 
is  because  of  the  limited  aperture  of  the  receiver  array. 

Figures  6c  and  d  give  the  F-K  spectra  for  the  seismograms  recorded  over  the  2-D  and  3-D 
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interfaces.  The  F-K  spectra  for  the  2-D  interface  clearly  shows  that  similar  amounts  of  both 
backscattered  and  forward  scattered  energy  are  present  in  the  seismogram.  This  distribution  is 
in  agreement  with  the  numerical  mean  reflection  coefficients  presented  for  the  2-D  interface  by 
Schultz  and  Toksoz  (1993a),  which  show  similar  amounts  of  both  forward  and  back  scattering.  The 
structure  of  the  spectra  is  parallel  to  the  Ky  just  like  the  rough  surface  is  parallel  to  the  Y  axis. 
The  F-K  spectra  for  the  3-D  interface  also  shows  a  large  amount  of  both  forward  and  back  scattered 
energy.  However,  the  forward  scattered  energy  appears  to  be  dominant.  This  suggests  that  more 
forward  scattering  appears  to  exist.  However,  one  must  remember  that  it  is  possible  for  energy, 
which  is  diffracted  back  towards  the  seismic  line  from  a  side  scatterer,  to  appear  with  a  positive 
Kx  wavenumber.  Thus,  in  the  3-D  case,  out-of-plane  scattering  can  give  the  appearance  of  more 
forward  scattering  and  less  backscattering.  Nevertheless,  these  results  appear  to  support  the  results 
of  Schultz  and  Toksoz  (1993b),  which  showed,  based  on  experimental  data,  that  3-D  interfaces  tend 
to  show  a  more  rapid  decrease  in  backscattering  combined  with  a  more  rapid  increase  in  forward 
scattering  as  an  incident  wave  arrives  at  larger  incident  angles  (with  respect  to  the  vertical).  In 
the  3-D  case  the  wide  spread  of  energy  in  the  F-K  spectra  clearly  demonstrates  the  out-of-plane 
reflections. 

Figure  7  gives  the  rms  amplitude  of  the  primary  reflection  versus  incident  angle.  The  seismogram 
section  is  taken  from  Y=l,  where  the  3-D  and  the  2-D  surface  have  the  same  profile.  The  geometric 
spreading  is  compensated  according  to  the  flat  interface.  The  rms  amplitudes  of  the  first  cycle  is 
shown  in  Figure  Ta.The  2-D  and  3-D  rms  amplitudes  are  similar,  but  they  are  both  smaller  than 
the  flat  interface  amplitude.  This  is  exactly  what  we  expected.  As  we  add  more  cycles  into  the 
rms  amplitude  (Figure  7b),  more  scattered  energy  is  accounted  for.  It  clearly  shows  that  the  extra 
degree  of  freedom  associated  with  the  3-D  interface  geometry  can  result  in  a  larger  contribution  of 
scattered  energy.  The  3-D  variation  of  the  interface  topography  results  in  a  larger  amplitude  than 
the  2-D  one  near  normally  incident  angles.  This  gives  a  first  estimate  of  the  scaling  differences 
between  the  2-D  and  3-D  reflection  coefficients  studied  by  Schultz  (1994)  in  his  thesis,  where  the 
rms  amplitude  from  2-D  numerical  calculations  and  3-D  ultrasonic  experiments  are  compared.  The 
out-of-plane  contributions  from  3-D  topography  also  result  in  the  amplitude  being  larger  than  the 
flat  interface. 
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SUMMARY 


The  3-D  time  domain  finite  difference  method  is  used  to  investigate  the  rough  fiuid-solid  interface 
scattering.  The  scheme  is  implemented  on  the  nCUBE  2S  massively  parallel  computer.  The 
synthetic  reflection  data  is  analyzed  by  the  F-K  spectra.  The  2-D  scattering  clearly  shows  the 
similar  amount  of  forward  and  backward  scattering,  while  in  the  3-D  case  out-of-plane  reflections 
give  more  forward  scattering  and  less  back  scattering.  The  3-D  effects  are  also  shown  on  the  larger 
reflected  rms  amplitude  than  the  2-D  case  near  normally  incident  angles.  This  is  because  of  the 
contributions  from  out-of-plane  scattering  from  3-D  topography.  This  3-D  out-of-plane  scattering 
can  be  easily  identified  on  the  F-K  spectra. 

In  this  paper  we  only  take  a  first  step  to  demonstrate  the  difference  of  2-D  and  3-D  surface 
scattering.  More  thorough  investigations  are  needed.  For  example,  we  need  to  test  a  large  amount  of 
realizations  of  the  interface  to  obtain  statistical  insight  into  2-D  and  3-D  rough  surface  scattering. 
In  our  simulation  we  only  consider  the  point  source.  To  achieve  high  angular  resolution  of  the 
reflection,  we  need  to  form  a  beam  to  the  interface  instead  of  the  point  source  (Stephen  and  Swift, 
1994).  This  will  give  us  a  more  quantitative  measure  of  the  2-D  and  3-D  reflection. 
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P  wave  velocity  a  S  wave  velocity  /?  density  p 
(m/s)  (m/s)  (g/c.c.) 

Fluid  1500  —  1.0 

Solid  4000  2300  2.3 

Table  1:  The  velocity  and  the  density  values  of  the  fluid  and  the  solid. 
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Figure  1:  Staggered  grid  used  to  discretize  equations  (2.5)  and  (2.6).  Solid  circles 
represent  the  velocities.  Open  circles  represent  the  shear  stresses.  The  solid  square 
represents  the  normal  stresses. 
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Figure  2:  C 
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Figure  4:  One  example  of  a  Gaussian  rough  surface.  The  correlation  length  is  100  m 
and  the  rms  slope  is  30  degrees,  (a)  2-D  rough  interface,  (b)  3-D  rough  interface. 
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(b)  3-D  Rough  Surface 
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Figure  4:  continued 
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(a)  Flat  Interface 
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Figure  5:  Synthetic  seismograms  of  fluid-solid  interface  reflections.  They  are  2-D  slices 
of  the  3-D  data  set.  The  direct  arrivals  are  removed  from  the  data  set.  (a)  flat  interface, 
(b)  2-D  rough  interface,  (c)  3-D  rough  interface. 
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Figure  5:  continued 
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(c)  3D  Interface 
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Figure  5:  continued 
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(a)  Direct  Wave 


Figure  6:  FK  plots  of  synthetic  data  at  frequency  15  Hz.  (a)  direct  arrivals,  (b)  flat 
interface  reflection,  (c)  2-D  rough  interface  reflection,  (d)  3-D  rough  interface  reflection. 
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Figure  6:  continued 
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(c)  2D  Interface 
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Figure  6:  continued 
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(d)  3D  Interfac 
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Figure  6:  contiu’ 


Source  time  functions  of  nuclear  explosions  and  earthquakes  in 
Central  Asia  determined  using  empirical  Green’s  functions 

Yingping  Li,  M.  Nafi  Toksoz,  and  William  Rodi 
Earth  Resources  Laboratory 

Department  of  Earth,  Atmospheric,  and  Planetary  Sciences 
^Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139,  USA 

ABSTRACT 

Relative  source  time  functions  (RSTF)  have  been  estimated  for  four  underground 
nuclear  explosions  and  seven  earthquakes  in  Central  Asia  using  broadband  P 
waveforms  of  nearby  smaller  events  as  empirical  Green’s  functions  (EGF).  RSTFs 
of  the  four  explosions  {mb  =  5.3  to  6.5)  are  each  dominated  by  a  simple  pulse  with 
a  source  duration  of  0.4  to  0.8  s.  RSTFs  for  two  of  the  explosions  show  a 
significant  secondary  pulse  with  a  pulse  width  similar  to  that  of  the  first  pulse.  We 
conclude  that  the  secondary  phases  are  most  likely  associated  with  the  spall 
slapdown  phenomenon.  Seismic  moment  releases  by  the  spall  phases  are  less  than 
one  third  of  those  by  the  first  explosion  pulses.  Elastic  radii  of  the  explosions  are 
estimated  to  be  about  0.25  to  0.5  km,  and  stress  drops  of  the  explosions  range  from 
13  to  52  MPa.  In  contrast,  RSTFs  of  earthquakes  studied  {mb  =  5.5  to  6.6) 
typically  comprise  multiple  source  pulses  with  a  total  source  duration  from  a  few  to 
several  tens  of  seconds,  indicating  that  the  complex  source  process  involves  a  fault 
dimension  of  several  tens  of  kilometers.  Stress  drops  of  the  earthquakes  are  much 
smaller  than  those  of  the  nuclear  explosions,  ranging  from  0.5  to  3.5  MPa.  Our 
study  demonstrates  the  power  of  the  EGF  method  for  retrieving  RSTFs  and  reveals 
that  differences  in  RSTFs  and  source  parameters  can  be  used  to  distinguish  large 
nuclear  explosions  from  moderate  to  large  earthquakes  {mb  ^  5.5). 
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INTRODUCTION 


The  seismic  source  time  function  (STF)  of  an  underground  nuclear  explosion  contains 
important  parameters  that  characterize  the  complex  physical  processes  that  occur  during  the 
explosion  [e.g.,  Masse,  1981;  Denny  and  Johnson,  1991].  During  the  last  three  decades, 
seismologists  have  developed  different  methods  to  characterize  the  STFs  and  other  source 
parameters  of  nuclear  explosions  in  the  time  or  the  frequency  domain  and  to  compare  them 


with  those  of  natural  earthquakes  [e.g.,  Brune  and  Pomeroy,  1963;  Toksdz  et  al.,  1964, 
\965',  Haskell,  1961',  Mueller  and  Murphy,  \91\',  Helmberger  and  Harkrider,  1972;  von 
Seggem  and  Blandford,  1972;  Muller,  1973;  et^al,  1914;  Murphy,  1977;  Burdick  and 


Helmberger,  1979;  Helmberger  and  Hadley,  1981;  iMy  et  ai,  1984a;' Der  et  al.,  1987; 
Denny  and  Goodman,  1990;  Douglas,  1991;  Patton,  1991;  Stump  and  Reinke,  1991]. 

Inference  of  the  STF  of  a  nuclear  explosion  from  observed  seismograms  requires 
isolating  the  source  effects  from  the  effects  of  wave  propagation  through  the  Earth, 
recording  site  response,  and  instrument  response.  Burdick  and  Helmberger  [1979] 
determined  effective  STFs  for  several  megaton-class  nuclear  explosions  by  deconvolving 
the  theoretical  instrument  response  and  Q  operator.  Der  et  al.  [1987]  developed  a  multiple- 
channel  iterative  deconvolution  procedure  to  obtain  effective  STFs  by  removing  a  receiver 
structure  response  and  the  instrument  response.  Stump  and  Reinke  [1991]  used  theoretical 
Green’s  functions  to  account  for  propagation  effects  and  carried  out  a  moment  tensor 
inversion  to  determine  the  STF  for  explosion  COALORA  (<  20  kt).  The  success  of  these 
methods  depends  largely  upon  detailed  knowledge  of  the  Earth's  structure  and  attenuation 
effects.  As  an  alternative,  we  shall  remove  the  responses  of  the  Earth’s  structure  and  the 
seismic  instrument  from  observed  seismograms  and  estimate  the  STFs  for  several 
explosions  in  Central  Asia  by  using  the  empirical  Green’s  function  (EGF)  method.  The 
great  advantage  of  this  approach  is  that  detailed  knowledge  about  the  Earth's  structure  and 
attenuation  is  not  required. 

The  essential  idea  behind  the  EGF  approach  is  as  follows.  For  two  seismic  events 
having  a  similar  hypocenter  and  focal  mechanism  but  different  sizes  (e.g.,  two  nuclear 
explosions  at  the  same  test  site,  or  an  earthquake  and  its  foreshock  or  aftershock),  one  can 
treat  the  waveform  of  the  smaller  event  as  EGF  [Hartzell,  1978]  and  deconvolve  it  from 
that  of  the  larger  event  to  obtain  a  relative  STF  [Mueller,  1985].  Since  the  pair  of  events 
shares  the  same  propagation  path  and  is  recorded  at  the  same  station,  the  deconvolution 
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procedure  should  result  in  a  relative  STF  that  is  corrected  for  all  path,  site,  and  instrument 
effects.  The  retrieved  STF  represents  the  far-field  displacement  pulse,  which  is 
proportional  to  the  time  derivative  of  seismic  moment  {Aki  and  Richards,  1980].  Therefore 
the  STF  we  use  here  is  the  reduced  velocity  potential  (RVP)  rather  than  the  reduced 
displacement  potential  (RDP),  which  is  commonly  called  the  seismic  source  function  in  the 
literature  of  explosion  sources  [De/iny  amf  /oftn.son,  199 1]‘. 

The  EGF  method  has  been  extensively  applied  to  estimate  STFs  and  to  retrieve  source 
parameters  for  microearthquakes  and  large  earthquakes,  using  local  seismic  network  data 
[t.g.,  Mueller,  \9%5‘,  Frankeh  et  al.;  19S6‘,  Hutchings,  1987;  Li  and  Thurber,  1988; 
Hutchings  and  Wu,  1990;  Mori  and  Frankel,  1990;  Mori  and  Hartzell,  1990;  Xie  et  al, 
1991;  Hough  et  al.,  1991;  Mori,  1993],  strong  motion  data  [e.g.,  Hartzell,  1978;  Chen  et 
al,  1991;  Kanamori  etal,  1992],  and  regional  and  teleseismic  waveforms  [e.g.,  Weidner 
and  Aki,  1973;  Patton,  1980;  Hartzell,  1989;  Ammon  et  al,  1993;  Li  and  Toksdz,  1993; 
Velasco  et  al,  1994].  Burger  and  Langston  [1985]  also  applied  the  EGF  method  to  study 
the  source  mechanism  of  the  May  18,  1980,  Mount  St.  Helens  eruption  from  regional 
surface  waves.  However,  few  endeavors  had  been  made  to  extract  the  STFs  of  nuclear 
explosions  with  the  EGF  deconvolution  method  until  recently  [Li  et  al,  1993;  Toksdz  et 
al,  1993;  Jiao  and  Wallace,  1993].  In  this  paper,  we  shall  retrieve  the  relative  STFs  of 
four  nuclear  explosions  and  seven  earthquakes  in  Central  Asia  (Figure  1)  using  the 
teleseismic  and  regional  broadband  P  waveforms  of  the  smaller  events  with  similar 
locations  as  the  EGFs  and  compare  the  RSTFs  of  explosions  with  those  of  natural 
earthquakes  to  examine  the  potential  value  of  this  method  in  nuclear  monitoring  and 
discrimination. 


DATA  AND  METHODS 

The  locations  of  nuclear  explosion  pairs  in  the  Kazakhstan  Test  Site  (KST)  at  Balapan, 
Kazakhstan,  and  in  the  Xinjiang  Test  Site  (XTS)  at  Lop  Nor,  Xinjiang,  China,  and  seven 
earthquake  pairs  in  Central  Asia,  as  well  as  locations  of  some  of  the  seismic  stations  used 
in  this  study  are  shown  in  Figure  1.  Tables  1  and  2  present  hypocentral  parameters  and 
magnitudes  from  the  U.S.  Geological  Survey  (USGS)  Preliminary  Determination  of 
Epicenters  (PDE)  for  the  explosions  and,  earthquakes.  The  depths  of  burial  (DOB)  for 
explosions  in  Table  1  were  estimated  with  the  formula  of  Jih  and  Wagner  [1991], 
log(DOB)  =  0.31/71/7  +  0.835.  The  vertical  component,  very  broadband  (VBB)  waveforms 
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used  in  this  research  were  collected  by  the  Incorporated  Research  Institutions  for 
Seismology,  Global  Seismographic  Network  (IRIS  GSN),  Joint  Seismic  Program  (ISP) 
and  Chinese  Digital  Seismographic  Network  (CDSN)  and  were  obtained  through  the  IRIS 
Data  Management  Center  (DMC).  The  IRIS  instrument  system  amplitude  response  is 
typically  flat  to  velocity  from  360  s  to  10  Hz  and  the  sampling  rate  is  20  samples  per 
second  for  the  VBBchanneL'’  I  ■  .i  i  :  i 

Application  of  the  EGF  method  requires  the  hypocentral  separation  between  the  larger 
event  and  the  EGF  event  to  be  much  smaller  than  the  distance  to  a  given  receiver  and 
comparable  to,  or  less  than,  the  shortest  wavelength  being  analyzed,  so  that  the  two  events 
share  the  same  wave  path  to  a  given  receiver  except  for  small  differences  in  the  source 
region  [e.g.,  Patton,  1980],  Based  on  the  PDE  locations  for  two  pairs  of  explosions  at 
KTS  (Figure  2a),  the  separation  between  the  881112  (read  as  November  12,  1988;  year, 
month,  day)  explosion  {mb  =  5.3)  and  the  880614  explosion  (EGF  event,  tnb  =  5.0)  is 
about  4  km;  while  the  distance  between  the  880914  Joint  Verification  Experiment  (JVE) 
explosion  {mb  =  6.1)  and  the  890708  explosion  (EGF  event,  mb  =  5.6)  is  about  7  km. 
However,  SPOT  satellite  image  analysis  [Thurber  et  ah,  1993]  indicates  that  the 
separations  between  these  event  pairs  are  approximately  1.5  and  3.5  km,  respectively 
(Figure  2b).  PDE  locations  for  three  explosions  at  XTS  are  shown  in  Figure  2c.  The 
distance  between  the  920521  explosion  {mb  =  6.5)  and  the  900816  explosion  (EGF  event, 
mb  =  6.2)  is  less  than  5  km,  and  the  separation  between  the  mb  =  6.5  event  and  the  931005 
explosion  (EGF  event,  mb  =  5.9)  is  estimated  to  be  12  km.  Unfortunately,  no  satellite 
image  analysis  is  available  for  these  events.  To  estimate  the  location  separation 
independently,  we  fixed  the  explosion  depths  at  the  surface  and  used  differences  between 
R-P  times  (lags  between  P  and  Rayleigh  waves)  among  the  events  from  a  few  stations  to 
derive  the  locations  relative  to  the  largest  explosion.  The  cross-correlation  technique  [e.g., 
Pechmann  and  Kanamori,  1982;  Poupinet  et  al,  1984]  is  applied  to  P  and  Rayleigh  wave 
windows  for  obtaining  precise  readings  of  relative  arrival  times  of  the  waves  from  the  VBB 
seismograms.  Reading  error  is  typically  less  than  ±0.05  s,  resulting  in  a  location  error  less 
than  0.5  km.  The  relocations  (Figure  2d)  indicate  that  the  distance  between  the  920521  and 
900816  explosions  is  about  2  km,  and  the  distance  between  the  920521  and  931005 
explosions  is  about  7  km. 

Assuming  two  nearby  surface  sources  with  epicenter  separation  AL  in  the  N-S 
direction,  the  difference  between  two  ray  path  lengths  (AD)  for  the  two  events  can  be 
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simply  written  as  AD  =  zyL(sin(|)cos0),  where  ^  is  the  takeoff  angle  and  q  is  the  source-to- 
receiver  azimuth.  For  an  event  pair  recorded  at  a  distance  of  35°  and  9  =  0°,  the  epicentral 
separation  of  7  km  will  correspond  to  the  maximum  difference  between  two  path  lengths  of 
about  2.8  km.  Given  an  average  P  wave  velocity  in  KTS  and  XTS  of  6.5  km/s  [Quin  and 
Thurber,  1992;  Mangino  and  Ebel,  1992]  and  the  dominant  frequencies  of  the  explosion 
STFs  ranging  from  0.3  to  3  Hz,  the  shortest  wavelength  we  fesdlVe  will  be  about  2.3  km, 
which  is  comparable  to  AD.  For  teleseismic  and  regional  data  used  here,  we  find  that 
events  with  source  separations  up  to  7  km  are  still  similar  to  each  other  and  can  be  used  as 
EGFevents.^  '  '  ' 

To  illustrate  the  EGF  deconvolution  method,  figure  3  shows  vertical  component 
broadband  P  waveforms  of  two  XTS  explosions  (920521  and  900816)  recorded  at  station 
OBN,  about  4000  km  from  the  test  site.  We  note  that  the  first  wavelets  of  the  two 
explosions  are  very  similar,  reflecting  intercorrelation  [Lay  et  ai,  1984b]  between  the 
nearby  explosions.  The  observed  seismogram  Ul(t)  of  the  larger  explosion  can  be 
expressed  as  the  convolution  of  its  STF,  5/(0,  the  impulse  response  of  path  P{t),  the 
receiver  structure  R(t),  and  the  instrament  response  I(t): 


Ui(t)  =  Slit)  *  Pit)  *  Rit)  *  1(0, 


(1) 


where  asterisk  represents  convolution.  For  an  EGF  event  with  its  STF,  5g(0,  its 
seismogram  Ug(t)  can  be  written  as 

Ugit)  =  Sgit)  *  Pit)  *  Rit)  *  7(0.  (2) 


The  broadband  P  waveform  of  the  smaller  event  is  used  as  the  EGF  and  deconvolved  from 
that  of  the  larger  explosion  to  obtain  a  relative  STF  for  the  larger  event.  The  deconvolution 
in  the  frequency  domain  is  a  spectral  division: 


f/g(CO)  5g(C0) 


(3) 


where  (O  is  the  angular  frequency  and  Sr{cd)  is  the  spectrum  of  the  RSTF  of  the  larger 
event.  The  deconvolution  procedure  used  here  is  similar  to  that  of  Li  and  Thurber  [1988] 
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and  Li  and  Toksdz  [1993].  The  P  waveforms  of  both  events  are  windowed  to  20  s,  zero 
padded  to  8192  points,  and  then  Fourier  transformed.  The  spectrum  of  the  large  explosion 
is  divided  by  that  of  the  EGF  event  and  then  inverse  transformed  to  obtain  the  RSTF,  Sr(t). 
To  stabilize  the  deconvolution,  we  used  a  Gaussian  filter  to  slightly  smooth  the  spectra 
before  the  spectral  division.  After  deconvolution,  a  fourth-order  Butterworth  low-pass 
filter,  with  a  comer  frequency  of  3  Hz;  iis  used  td'reduce  the  high-frequency  noise.  The 
RSTF  of  the  920521  explosion  comprises  a^sirnple  pulse  followed  by  a  secondary  phase, 
with  a  total  source  duration  of  about  1.6  s  (Figure  3b).  We  look  back  at  Figure  3a  and  find 
a  secondary  arrival  on  the  seismog'rkin  of  the  lafget  event:  'The  secondary  phase  is  possibly 
caused  by  spall  slapdown  or  tectonic  release.  We  will  discuss  this  matter  in  detail  later. 

The  deconvolved  pulse  is  the  RSTF,  and  its  pulse  width  presents  the  relative  source 
duration  of  the  larger  event  compared  to  that  of  the  EGF  event  [Frankel  et  ai,  1986].  The 
tme  STF  of  the  larger  event,  Si(t),  will  be  the  convolution  of  the  source  pulse  of  the  EGF 
event  with  the  RSTF  of  the  larger  event 

Si(t)^Sr(t)*Sg(t).  (4) 

Si(t)  =  Srit)  only  when  Sg{t)  =  d(r).  If  the  source  duration  of  the  EGF  event  is  short 
enough  to  approximate  a  delta  function,  then  the  RSTF  can  be  approximately  equal  to  the 
true  STF.  However,  if  the  source  duration  of  the  EGF  event  is  nonnegligible,  larger  error 
will  be  introduced  for  estimating  the  source  duration  from  the  RSTF,  especially  for  the 
events  which  have  relatively  narrow  deconvolved  pulse  width  compared  to  the  STF  pulse 
width  of  the  EGF  event.  We  can  obtain  the  true  STF  of  the  larger  event  from  equation  (4) 
if  the  STF  of  EGF  event  is  known. 

RELATIVE  SOURCE  TIME  FUNCTIONS  FOR  EXPLOSIONS 

With  broadband  P  waveforms  recorded  at  station  WMQ  (A  =  9°),  we  estimated  the  RSTF 
of  the  881112  explosion  (m^  =  5.3)  in  KTS  using  the  880614  explosion  (mb  =  5.0)  as  the 
EGF  event  (Figure  4a).  The  RSTF  of  the  881112  explosion  is  a  simple  pulse  with  a 
source  duration  of  0.4  s,  suggesting  a  rapid  energy  release  of  the  nuclear  explosion.  We 
attributed  the  success  of  retrieving  the  RSTF  to  the  selection  of  an  appropriate  EGF:  the 
spatial  separation  between  the  two  explosions  is  less  than  1.5  km  and  the  depth  difference 
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is  estimated  at  0.06  km  (Table  1  and  Figure  2).  Since  only  one  regional-distance  station 
was  used  to  retrieve  the  RSTF,  the  background  noise  level  in  Figure  4a  is  relatively  high. 
No  obvious  secondary  phase  can  be  resolved  from  the  single-station  estimate  of  the  RSTF. 

Broadband  P  waveforms  from  stations  WMQ,  LZH,  HIA,  and  COL  were  used  to  extract 
the  RSTF  of  the  880914  JVE  explosion  {mb  =  6.1)  using  the  890708  explosion  (mb  =  5.6) 
as  the  EGF  event.  The  results  are  shown  in  Figure  4b.  The  stacked  estimate  of  the  RSTF 
(bottom  trace  in  Figure  4b)  shows  a  major  pulse  with  a  total  duration  of  about  1.2  s.  The 
asymmetric  pulse  shape  of  the  RSTF  and  its  positive-polarity  shoulder  suggests  that  a 
secondary  pulse  with  a  much  smaller  amplitude  and  a  similar  pulse  width  follows  the  first 
explosion  pulse. 

The  RSTF  of  900816  explosion  (mb  =  6.2)  at  XTS  was  extracted  using  the  931005 
explosion  (mb  =  5.9)  as  EGF  event.  We  used  broadband  P  waveforms  recorded  at  stations 
OBN,  KEV,  and  COL  with  epicentral  distances  ranging  from  36°  to  65°.  The  single¬ 
station  estimates  and  the  stacked  STF  (bottom  trace)  are  presented  in  Figure  5a.  The  major 
feature  of  the  STF  is  a  simple  pulse  with  a  duration  of  0.55  s.  No  significant  secondary 
phase  can  be  clearly  identified  on  the  stacked  trace.  However,  we  do  observe  the  variation 
of  the  STFs  among  the  single-station  estimates,  especially  for  the  latter  part  of  the  STFs. 
This  obser\'ation  can  be  interpreted  as  either  the  tectonic  release  effect  or  the  imperfection  of 
the  EGF  caused  by  a  relatively  larger  epicentral  separation  between  the  two  events. 

We  retrieve  the  RSTF  of  the  920521  explosion  (mb  =  6.5)  at  XTS  with  the  broadband  P 
waveforms  of  two  nearby  explosions  (931005,  mb  =  5.9  and  900816,  mb  =  6.2)  as  the 
EGFs.  This  data  set  provides  an  opportunity  not  only  to  extract  the  RSTF  with  more  than 
one  EGF  event  but  also  to  examine  the  conditions  under  which  the  EGF  method  can  be 
successfully  applied.  The  RSTF  extracted  using  the  931005  explosion  (mb  =  5.9)  as  the 
EGF  event  is  shown  in  Figure  5b,  while  the  RSTF  retrieved  using  the  900816  explosion 
(mb  =  6.2)  as  the  EGF  event  is  presented  in  Figure  5c.  The  RSTFs  of  the  920521 
explosion  are  dominated  by  a  simple  pulse  with  a  duration  about  0.8  s,  followed  by  a 
significant  secondary  pulse  with  a  smaller  amplitude  and  a  similar  pulse  width.  The  total 
source  duration  of  the  explosion  is  about  1.6  s.  The  area  under  the  RSTF  is  proportional  to 
the  seismic  moment  release  [Aki  and  Richards,  1980],  and  the  seismic  moment  release  of 
the  mb  =  6.5  explosion  is  larger  than  those  of  two  EGF  event  (931005  and  900816)  by 
factors  of  5.5  and  2.5,  respectively.  This  result  agrees  with  the  fact  that  the  seismic 
moment  release  ratio  between  the  900816  and  931005  explosion  is  about  2.2,  which  is 
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estimated  using  the  area  under  the  RSTF  in  Figure  5a,  indicating  that  the  EGF  method  can 
provide  a  reliable  means  to  estimate  the  relative  seismic  moment  release. 

In  spite  of  having  used  two  events  with  different  sizes  {mb  =  5.9  and  6.2)  and  different 
hypocenters  (depth  =  0.46  and  0.57  km,  AL  =  7  and  2  km)  as  EGF  events,  the  waveforms 
of  the  two  stacked  RSTFs  for  the  920521  explosion  are  almost  identical  except  for  the 
pulse  amplitudes  (Figures  5b  and  5c),  indicating  that  the  EGF  method  can  effectively 
correct  the  wave  propagation  effects  through  the  Earth  as  well  as  the  instmment  effect,  as 
long  as  one  or  more  suitable  EGF  events  are  available.  iThe  successful  extraction  of  the 
RSTF  using  the  two  EGF  events,  with  spatial  separations  of  about  2  and  7  km  to  the 
largest  explosion  {mb  =  6.5),  also  supports  our  earlier  argument  about  the  criteria  of 
selecting  suitable  EGF  events  for  applying  the  EGF  method  to  the  teleseismic  and  far- 
regional  P  waveform  data. 

With  the  mb  =  5.9  explosion  as  the  EGF  event,  inferred  RSTF  pulse  widths  are  0.55  s 
and  1.6  s  for  the  mb  =  6.2  and  mb  =  6.5  explosions,  respectively  (Figures  5a  and  5b). 
Thus,  for  explosions  with  a  size  difference  of  0.3  magnitude  units,  the  RSTF  pulse  width 
of  the  smaller  event  is  narrower  than  that  of  the  larger  event  by  a  factor  of  3.  If  this  is  a 
linear  relation  and  valid  for  the  events  with  smaller  magnitude,  the  RSTF  pulse  width  for 
the  =  5.9  explosion  is  expected  to  be  about  0.2  s.  A  scaling  relation  for  the 
earthquakes  in  continental  interiors  [Somerville  et  al,  1987]  indicates  that  the  source 
durations  are  expected  to  be  0.4  to  1.0  s  for  events  with  seismic  moments  of  5xl0l5  to 
5x1016  Nm.  The  source  duration  estimated  with  the  scaling  relation  is  comparable  to  our 
extrapolation  result,  although  the  scaling  relation  for  earthquake  may  not  be  entirely 
appropriate  for  the  explosions.  We  have  used  two  EGF  events  {mb  =  5.9  and  6.2)  to 
retrieve  the  RSTF  of  the  mb  =  6.5  explosion.  The  waveform  shapes  of  the  two  RSTFs  in 
Figures  5b  and  5c  are  very  similar,  and  the  pulse  width  difference  between  the  two  RSTFs 
is  hardly  observed.  This  result  implies  that  the  STF  pulses  of  two  EGF  events  {mb  =  5.9 
to  6.2)  may  have  relatively  short  source  duration  comparing  to  the  STF  of  the  largest 
explosion  and  that  the  source  pulse  of  an  EGF  event  with  a  duration  of  0.55  s  may  still  be 
approximately  treated  as  a  delta  function  in  this  case. 

To  examine  the  effects  of  source  pulses  of  EGF  events  on  the  true  STF  of  the  largest 
explosion,  we  convolve  five  assumed  triangular  source  pulses  of  the  EGF  event,  Sg{t), 
with  the  RSTF,  Sr{t),  in  Figure  5c  to  predict  the  true  STF,  5/(r),  for  the  largest  explosion 
(Figure  6a).  The  source  duration  of  Sg{t)  is  from  0.2  to  1.0  s  with  an  increment  of  0.2  s 
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and  the  peak  amplitudes  of  the  triangular  pulses  are  10,  5,  0.33,  2.5,  and  0.5,  so  that  the 
areas  under  the  STFs  are  one  unit.  The  convolved  STFs  show  a  little  variation  of  the  pulse 
width  of  the  predicted  STFs  (Figure  6a).  For  a  source  pulse  of  an  EGF  event  with  a 
duration  of  up  to  0.8  s,  the  pulse  width  of  the  predicted  STF  is  broadened  only  by  about 
0. 1  s.  This  result  indicates  that  the  RSTF  pulse  width  of  the  larger  event  is  relatively 
insensitive  to  the  pulse  width  variation  of  the  EGF,  events.  Therefore  the  RSTF  of  the 
larger  event  can  be  approximately  treated  as  the  STF  with  a  tolerable  error.  Figure  6b  plots 
the  source  spectra  for  the  RSTF,  i'rCco),  and  the  five  predicted  STFs.  The  corner 
frequencies,  /c,  of  the  five  predicted  STFs ,  range  froip  1.2  to  2  Hz  for  the  mb  =  6.5 
explosion. 

SOURCE  PARAMETERS  OF  EXPLOSIONS 

Figure  7  summarizes  the  RSTFs  for  the  four  explosions  in  Central  Asia.  With  the 
estimated  STFs  and  magnitudes  of  the  exploslibns,  one  can  derive  other  explosion  source 
parameters,  such  as  the  reduced  displacement  potential  (\|/„),  the  seismic  moment  {Mq),  the 

elastic  radius  (Re),  and  the  stress  drop  (As).  The  seismic  moment  of  an  explosion  is  given 
by  [Muller,  1973] 

Mo  =  4npaVoo.  (5) 

where  r  and  a  are  the  density  and  the  P  wave  velocity  at  the  source,  respectively.  Ringdal  et 
al.  [1992]  present  an  empirical  relationship  between  \|/„  and  mb  for  explosions  in  KTS, 

logvi/„=  l.lmt-2.57.  (6) 

Assuming  that  the  relation  is  also  valid  for  XTS  and  assuming  p  =  2700  kg/m^  and  a  = 
5.2  km/s  [Quin  and  Thurber,  1992;  Matzko,  1992],  we  can  calculate  reduced  displacement 
potentials  and  seismic  moments  for  all  explosions  (mb  =  5.0  to  6.5)  used  in  this  study 
(Table  3).  Since  the  area  under  the  RSTF  is  the  seismic  moment  ratio  (Mo/Mog)  between 
the  larger  explosion  and  the  EGF  event,  we  can  also  calculate  the  seismic  moment  of  the 
larger  explosion  by  multiplying  the  area  with  the  seismic  moment  of  EGF  event,  Mog-  In 
this  manner,  we  estimated  the  seismic  moment  not  only  for  the  explosion  source  pulses  but 
also  for  the  secondary  phases.  The  seismic  moments  of  the  larger  explosions  derived  with 
the  two  methods  agree  with  each  other  quite  well  (Table  3). 
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The  elastic  radius  Re  is  estimated  with  the  relationship  derived  by  [Murphy,  1977] 

=  aiinf^,  (7) 

where /c  is  the  comer  frequency  of  the  STF  pulse  and  the  P  wave  velocity  at  the  source  is 
assumed  to  be,  5.2  km/s.  Assuming  the  source  pulse  width  of  the  EGF  event  is  less  than  a 
•half  of  pulse  width  of  the  RSTF,  we  can  estimate  the  comer  ifrequency,/c  ,^ and  its  error 
bar  from  the  RSTF  and  corrected  STFs  with  the  method  showing  in  Figure  6b.  The  stress 
drops  for  the  explosions  are  calculated  from  Aki  et  a/.’s  {1969]  formula  for  a  spherical 
explosion  source  '  ;  .■  "  i  i  ,  ■  < 

Aa=2iMjm%R].  (8) 

The  estimated  source  parameters  and  their  error  bars  are  summarized  in  Table  3.  The 
seismic  moments  range  from  1.7  x  10^5  to  3.3  x  10^6  Nm,  and  elastic  radii  range  from 
250  to  520  m  for  the  four  explosions  with  mb  =  5.3  to  6.5.  The  mean  values  of  stress 
drops  are  estimated  about  13  MPa  for  one  event  with  mb  =  5.3  and  about  27  to  52  MPa  for 
larger  explosions  (mb  >  6.0).  The  lower  and  upper  bounds  of  the  stress  drop  estimates  are 
5  and  191  MPa.  Aki  et  al.  [1969]  estimated  the  source  parameters  for  the  681219 
explosion  BENHAM  (mb  =  6.3)  in  NTS  and  found  a  seismic  moment  of  2.5x10^^  Nm 
and  elastic  radius  of  450  to  900  m,  resulting  in  an  estimate  of  the  stress  drop  of  40  to  300 
MPa.  Our  estimates  of  the  elastic  radii  and  stress  drops  for  explosions  at  KTS  and  XTS 
are  comparable  to  those  estimated  for  the  BENHAM  explosion. 

Spall  Slapdown  Phases 

A  significant  feature  of  the  RSTF  for  the  920521  explosion  is  that  it  comprises  two 
distinct  pulses  with  a  total  source  duration  of  1.6  s.  This  result  is  consistent  with  the  result 
for  the  same  explosion  extracted  by  Lay  and  his  colleague  (T.  Lay,  personal 
communication,  1993)  using  the  surface  wave  and  EGF  method.  The  cause  of  this 
secondary  pulse  is  an  interesting  problem.  Previous  studies  of  explosion  sources  have 
revealed  that  besides  the  explosion  process  itself,  two  other  physical  processes  contribute 
to  the  nonisotropic  components  of  an  explosion  STF:  spall  [e.g.,  Eisler  et  al,  1966; 
Frasier,  1972;  Viecelli,  1973;  Day  et  al,  1983;  Stump,  1985;  Schlittenhardt,  1991]  and 
tectonic  release  [e.g.,  Toksdz  and  Kehrer,  1972;  Wallace,  1991;  Patton,  1991].  A  pP 
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reflection  [e.g..  Lay,  1991]  may  also  affect  the  RSTF,  since  significant  depth  difference 
between  an  explosion  and  its  EGF  event  may  cause  a  secondary  pulse  in  the  RSTF.  We 
will  try  to  interpret  the  secondary  pulse  from  all  possible  aspects. 

Assuming  that  the  DOB-m^  relation  [Jih  and  Wagner,  1991]  is  also  valid  for 
explosions  at  XTS,  then  the  depth  differences  between  the  mb  =  6.5  explosion  (920521) 
and  the  two  EGF  events  (931005  and  900816)  estimated  using  mb  (Table  1)  are  only  0.25 
and  0.14  km,  respectively.  If  the  secondary  phase  is  caused  by  the  depth  phase  pP,  for  an 
average  P  wave  velocity  of  5.2  km/s  for  the  shallow  crust  at  XTS  [Matzko,  1992],  the  time 
delay  of  0.75  s  between  the  two  pulse  peaks  will  correspond  to  a  depth  difference  of  about 
1 .95  km,  which  is  larger  than  the  ones  estimated  from  mb  by  a  factor  of  8  to  14.  The  delay 
time  of  0.75  s  is  too  large  to  associate  with  the  effect  of  the  depth  phase.  Thus  the 
secondary  pulse  was  probably  caused  by  spall  or  tectonic  strain  release,  rather  than  by  the 
pP  reflection.  For  the  same  reason,  we  believe  that  the  secondary  pulse  of  the  JVE 
explosion  (Figure  4b)  is  not  the  effect  of  the  pP  reflection  either. 

The  tectonic  release  caused  by  an  explosion  typically  generates  a  non-isotropic  seismic 
radiation  pattern  like  a  double-couple  earthquake  [e.g.,  Toksdz  and  Kehrer,  1972;  Wallace 
et  ai,  1985].  If  the  RSTF  of  the  explosion  (920521,  mb  =  6.5)  does  indeed  include  the 
component  of  tectonic  release,  it  is  expected  that  a  radiation  pattern  similar  to  an  earthquake 
will  be  observed.  We  looked  at  the  single  station  estimates  of  RSTFs  at  different  azimuths 
in  Figure  5c  and  found  only  slight  differences  among  the  RSTFs.  The  positive  polarities  of 
the  secondary  pulses  at  all  stations  (Figure  5c)  led  us  to  conclude  that  the  secondary  pulse 
is  dominated  by  the  contribution  from  the  spall  slapdown  phase,  which  is  thought  to  be 
associated  with  a  vertical  dipole  [e.g..  Stump  and  Reinke,  1991].  Since  the  tectonic  release 
component  of  an  explosion  can  be  best  observed  from  the  S  wave  and/or  surface  wave  data 
[Wallace,  1991],  we  suggest  that  a  further  study  with  S  wave  and  surface  waves  of  the 
same  explosions  will  lead  to  a  better  understanding  of  the  tectonic  release  components  of 
the  explosions. 

The  seismic  moment  release  of  the  secondary  pulse  from  the  880914  KST  and  the 
920521  XTS  explosions  is  about  10%  to  30%  of  that  for  the  first  explosion  pulse  (Table 
3).  It  is  interesting  to  compare  the  features  of  the  spall  slapdown  phases  we  observed  for 
the  explosions  at  XTS  and  KTS  with  those  of  the  spall  slapdown  phases  generated  by  other 
explosions  at  different  test  sites.  Table  4  summarizes  the  parameters  for  spall  slapdown 
phases  of  explosions  at  different  test  sites  estimated  with  different  methods  by  many 
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investigators.  For  explosions  {mb  =  3.8  to  6.3)  in  Nevada  Test  Site  (NTS),  the  delay 
times  between  the  spall  slapdown  phases  and  the  P  arrivals  {Pss-P)  range  from  0.17  to  1.7 
s,  and  the  pulse  widths  of  the  spall  phases  vary  from  0.3  to  1 .0  s.  For  the  large  explosions 
LONGSHOT,  MILROW,  and  CANNIKIN  {mb  =  6  to  7)  in  Amchitka  Test  Site  (ATS), 
Pss-P  ranges  from  0.7  to  2  s,  with  a  mean  of  1.35  s,  and  the  spall  pulse  widths  range 
from  0.5  to  I  s.  These  numbers  are  comparable  to  what  we  have  observed  for  the  spall 
slapdown  phase  generated  by  the  920521  explosion  and  the  880914  explosion. 

,  RELATIVE  SOURCE  TIME  FUNCTIONS  FOR  EARTHQUAKES 

The  RSTF  of  the  931002  earthquake  {mb  =  6.2)  in  southern  Xinjiang  was  estimated 
using  the  teleseismic  and  regional  broadband  P  waveforms  of  a  nearby  aftershock  {mb  = 
5.7)  as  the  EGFs  (Figure  8a).  Eight  single-station  estimates  of  the  RSTF  are  stacked  to 
obtain  an  average  RSTF.  The  mainshock  consists  of  two  subevents:  a  smaller  precursor 
event  followed  by  a  larger  one  2  s  later.  The  total  relative  source  duration  of  the  mainshock 
is  about  9  s.  The  location  of  the  earthquake  pair  is  about  300  km  south  of  XTS. 

We  retrieved  relative  RSTFs  for  a  deep  earthquake  (930809,  mb  =  6.3)  in  the  Hindu 
Kush  region  using  the  teleseismic  broadband  P  waveforms  of  its  foreshock  {mb  =  5.8)  as 
the  EGFs  (Figure  8b).  Two  pulses  can  be  clearly  observed  on  the  single  and  the  stacked 
RSTF  traces.  It  is  evident  that  the  mainshock  was  a  double  event  with  a  total  relative 
source  duration  of  about  14  s.  Both  mainshocks  (931002  and  930809)  show  a  relative 
complex  source  process  and  a  longer  relative  source  duration  compared  with  the 
explosions. 

We  also  apply  the  EGF  method  to  the  other  five  pairs  of  earthquakes  in  Central  Asia 
(Figure  1  and  Table  2).  Figure  9  summarizes  the  stacked  relative  STFs  estimated  for  all 
seven  mainshocks.  Stations  used  in  retrieving  RSTFs  for  each  event  are  listed  in  Table  5. 
The  body  wave  magnitudes  {mb)  for  the  mainshocks  range  from  5.5  to  6.6,  while  the 
surface  wave  magnitudes  {Ms)  are  from  5.2  to  7.4.  The  relative  source  duration  of  the 
earthquakes  is  from  7  to  29  s.  It  can  be  seen  clearly  from  the  RSTF  in  Figure  9  that  almost 
all  of  the  earthquakes  consist  of  multiple  events,  indicating  the  complicated  source 
processes  of  the  earthquakes  in  Central  Asia.  The  simplest  RSTF  among  the  seven 
earthquakes  is  900614  Kazakhstan  mb  =  6.1  earthquake.  Its  RSTF  contains  one  major  and 
one  minor  pulse  with  a  total  duration  of  about  7  s.  The  most  coniplicated  RSTF  is  for  the 
April  26,  1990,  Qinghai,  China,  mainshock  sequence.  The  PDE  reported  that  the  sequence 
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included  three  subevents  with  mb  =  5.7,  6.5,  and  6.2.  The  three  major  pulses 
corresponding  to  the  three  events  can  be  clearly  identified  from  the  RSTF  of  the 
mainshock. 

The  fault  radius  (r)  and  the  static  stress  drop  (Aa)  can  estimated  using  the  following 
relationships  1980;  Brwne,  1970]: 

,  - .  (9) 

1  -  (v^sin0/a) 

Aa  =  7Myl6r^  '  (10) 

where  x  1/2  is  the  rise  time  and  q  is  the  takeoff  angle;  vr  and  a  represent  rupture  velocity 
and  P  wave  velocity,  respectively.  We  use  a  =  6.5  and  8.2  km/s  and  v,-  =  2.0  and  2.5 
km/s  for  the  earthquakes  in  the  crust  and  the  upper  mantle,  respectively,  and  assume  an 
average  q  of  30°.  For  the  larger  event  in  each  pair,  the  seismic  moment  ratio  MolMog  is 
estimated  based  on  a  measurement  of  the  area  under  the  RSTF  pulse.  Using  the  CMT 
(HRV)  solutions  of  seismic  moments  for  the  EGF  events  (Mog)  as  reference,  we  can 
derive  the  seismic  moment  Mq  for  the  larger  events.  The  seismic  moments  Mo  estimated  in 
this  way  are  very  close  to  the  quick  CMT  solutions  [e.g.,  Ekstrom  et  al.,  1986]  for  the 
larger  events  (see  Table  6). 

As  we  mentioned  earlier,  the  traces  in  Figure  9  are  the  RSTFs  rather  than  the  STFs. 
The  RSTF  approximately  equals  the  STF  of  the  larger  event  only  when  the  source  pulse  of 
the  EGF  event  is  short  enough  to  approximate  a  delta  function.  For  the  EGF  events  with 
magnitudes  less  than  5.3,  the  RSTFs  could  be  approximately  treated  as  the  STFs  with 
tolerable  errors  since  the  source  pulse  widths  for  these  events  are  typically  1  to  2  s. 
However,  if  we  use  a  relatively  large  earthquake  (e.g.,  mb  >  5.7)  as  the  EGF  event  to 
retrieve  the  STF,  the  source  pulse  width  of  the  EGF  event  itself  must  be  accounted  for. 
Otherwise,  significant  errors  will  be  introduced  into  estimates  of  source  parameters.  We 
use  the  900614  Kazakh  earthquake  and  the  930809  Hindu  Kush  event  to  illustrate  the  effect 
of  source  pulse  width  of  EGF  event  on  the  RSTF  of  the  larger  event. 

According  to  the  quick  CMT  (HRV)  solution,  the  source  durations  for  the  900614 
Kazakh  mainshock  and  its  EGF  counterpart  (900803)  are  17  and  9  s,  respectively.  Our 
inferred  source  duration  of  the  RSTF  of  the  Kazakh  mainshock  is  only  7  s,  about  2.5  times 
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shorter  than  the  CMT  estimate.  We  assume  that  the  source  of  the  EGF  event  is  a  triangular 
pulse  with  duration  of  9  ±  1  s  and  convolve  them  with  the  RSTF  to  recover  an  "absolute" 
STF  (Figure  10a).  The  peak  amplitudes  of  these  triangle  pulses  are  0.25,  0.22,  and  0.2  so 
that  the  areas  under  the  pulses  are  one  unit.  The  duration  of  the  corrected  STF  is  measured 
about  18  ±  0.5  s,  close  to  the  CMT  estimate  for  the  larger  event.  Figure  10b  shows  how 
the  rise  time  xi/2  and  its  error  bars  are  measured  from, the  corrected  STFs  of  the  900614 
Kazakh  event.  The  source  durations  of  the  930809  Hindu  Kush  mainshock  and  its 
foreshock  are  16  and  7  s,  respectively,  based  on  the  CMT  (HRV)  solutions.  The  triangle 
source  pulses  with  duration  of  7  ±  1  s  was  convolved  with  the  RSTF  of  the  930809 
mainshock  to  retrieve  the  STF  of  the  mainshock  (Figure  10c).  The  source  duration  of  the 
corrected  STF  is  about  16.5  ±  1  s,  which  is  about  2  to  3  s  wider  than  that  of  the  RSTF. 
The  rise  time  and  measurement  error  for  the  first  subevent  of  the  930809  mainshock  are 
shown  in  Figure  lOd.  Using  the  same  approach,  we  measure  the  rise  times  for  the  931002 
and  920819  mainshocks  from  the  corrected  STFs.  For  the  900515,  900417,  and  900426 
mainshocks,  their  EGF  events  are  relatively  small,  from  magnitude  4.7  to  5.3.  We 
convolve  a  triangle  source  pulse  of  duration  1.5  ±  1  s  with  the  RSTFs  of  three  mainshocks 
and  find  that  the  corrected  STFs  differ  only  a  little  from  their  respective  RSTFs.  The  rise 
times  and  error  bars  estimated  from  the  "corrected"  STFs  for  12  subevents  of  the  seven 
mainshocks  are  listed  in  Table  6. 

Source  parameters  for  the  12  subevents  are  summarized  in  Table  6.  These  earthquakes 
have  seismic  moments  ranging  from  9.6  x  lO^^  to  7.7  x  10^^  Nm,  and  their  fault  radii  are 
estimated  at  about  4  to  21  km.  Estimated  stress  drop  values  fall  from  0.4  to  2.3  MPa.  The 
average  stress  drop  for  12  subevents  of  the  seven  mainshocks  is  1.4  ±  0.5  MPa.  It  should 
be  noted  that  the  stress  drop  estimates  are  made  under  the  assumption  of  a  rupture  velocity 
of  2  to  2.5  km/s.  If  the  rupture  velocity  is  assumed  to  be  3.0  to  3.6  km/s,  the  average 
stress  drop  will  decrease  to  0.5  ±  0.2  MPa.  However,  when  we  assume  a  relatively  low 
rupture  velocity  of  1.5  to  2  km/s,  the  average  stress  drop  will  increase  to  3.5  ±  1.3  MPa. 
Therefore  we  conclude  that  the  stress  drops  for  these  intraplate  earthquakes  in  Central  Asia 
are  typically  estimated  to  range  from  0.5  to  3.5  MPa. 

For  comparison,  we  also  estimate  the  rise  times  from  the  RSTFs  for  the  12  subevents. 
The  error  bars  for  these  relative  rise  time  measurements  are  typically  less  than  0. 1  s.  The 
source  parameters  estimated  with  the  relative  rise  times  are  also  listed  in  Table  6.  For  the 
RSTF  retrieved  using  relatively  small  events  {mh  <  5.3)  as  EGF  events,  stress  drops 
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estimated  from  the  RSTFs  differ  from  those  estimated  with  the  corrected  STFs  by  a  factor 
of  1.1  to  1.8,  suggesting  the  RSTFs  can  be  approximately  treated  as  the  STFs.  However, 
when  we  use  the  larger  earthquakes  (mb  =  5.7  to  6.0)  as  the  EGF  events  and  size 
differences  between  the  target  and  EGF  events  are  less  than  0.5  magnitude  units,  estimates 
of  the  rise  time  from  the  RSTFs  will  underestimate  the  fault  radius  by  a  factor  of  1.5  to  2.5 
and  overestimate  the  stress  drop  by  a  factor  of  3  to  16,  Therefore,  to  obtain  a  reliable 
estimate  of  source  parameters,  especially  stress  drop,  the  source  pulse  width  of  the  EGF 
event  should  be  taken  into  account  in  the  latter  case. 

The  lower  and  upper  bounds  of  the  stress  drop  estimates  for  the  earthquakes  in  Central 
Asia  are  0.4  and  21  MPa,  respectively.  The  wide  range  of  variations  indicates  the  stress 
drop  is  one  of  most  difficult  source  parameters  to  be  reliably  estimated.  The  error  sources 
include  the  uncertainties  in  the  assumed  rupture  velocity,  rise  time  measurement,  and 
circular  fault  assumption.  Our  preferred  stress  drop  estimates  of  0.5  to  3.5  MPa  are  much 
lower  than  those  for  some  moderate  earthquakes,  such  as  the  750204  mb  =  6.4  Haicheng, 
China,  earthquake  (14.4  MPa  [Chung  and  Brantley,  1989]),  the  831007  mb=  5.1 
Goodnow,  New  York,  earthquake  (26.5  MPa,  [Nabelek  and  Saurez,  1989]),  the  820109 
mb  =  5.8  New  Brunswick,  Canada,  earthquake  (50.0  MPa  [Nabelek,  1984]),  but 
comparable  to  those  of  other  moderate  earthquakes,  such  as  the  800727  mb  =  5.2 
Sharpsburg,  Kentucky,  earthquake  (0.3  to  0.6  MPa  [Herrmann  et  al,  1982]),  the  650912 
mb  =  5.9  South  China  Sea  earthquake  and  the  671210  =  6.0  Koyna,  India,  earthquake 

(1.8  and  2.0  MPa,  respectively  [Chung  and  Brantley,  1989]),  and  the  660628  mb  =  5.8 
Parkfield,  California,  earthquake  (3.3  MPa,  [Wu,  1968]). 

DISCUSSIONS  AND  CONCLUSIONS 

Analysis  of  the  RSTFs  of  four  nuclear  explosions  (mb  =  5.3  to  6.5)  and  seven 
earthquakes  (mb  =  5.5  to  6.6)  in  Central  Asia  indicates  that  the  earthquake  RSTF,  in 
general,  shows  a  complex  source  process  with  a  duration  of  several  tens  of  seconds, 
implying  that  the  earthquake  processes  involve  a  fault  dimension  of  a  few  to  several  tens  of 
kilometers.  In  contrast,  the  explosion  source  presents  a  relatively  simple  RSTF  with  one  or 
two  pulses  and  a  much  shorter  source  duration  of  about  0.4  to  1.6  s,  which  corresponds  to 
a  relatively  smaller  elastic  radius  of  0.25  to  0.5  km,  reflecting  the  rapid  energy  release  over 
a  small  volume  for  the  nuclear  explosions.  Among  the  four  explosions,  the  920521 
explosions  and  the  880914  explosions  show  significantly  a  secondary  pulse  with  a  pulse 
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width  similar  to  their  first  pulse.  We  interpret  that  the  secondary  pulses  are  most  likely 
associated  with  the  spall  slapdown.  The  moment  release  of  the  spall  slapdown  phase  is 
about  one  third  of  the  explosion  pulse  itself. 

The  source  durations  of  the  explosions  and  earthquakes  in  Central  Asia  are  plotted 
versus  the  magnitudes,  mb,  in  Figure  11.  For  events  with  similar  inagnitudes,  the  relative 
source  durations  for  explosions  are  shorter  than  those  of  e^hquake^  by  a  factor  of  about 
10,  indicating  the  earthquakes  and  explosions  within  this  magnitude  range  can  be  easily 
distinguished  from  one  another  by  examining  their  relative  source  durations.  Tables  3  and 
6  clearly  show  that  the  elastic  radii  of  the  explosions  are  much  smaller  than  the  fault 
dimension  of  earthquakes  with  similar  magnitudes  by  a  factor  of  10  to  40.  Stress  drops  of 
the  nuclear  explosion  sources  are  estimated  typically  at  13  to  52  MPa,  compared  with 
relatively  low  stress  drops  of  0.5  to  3.5  MPa  for  the  natural  earthquakes  with  similar 
magnitudes  in  the  vicinity  of  the  two  test  sites.  Our  results  strongly  indicate  that  the  large 
explosions  and  earthquakes  (mb  =  5.3  to  6.6)  in  Central  Asia  can  be  easily  distinguished 
from  each  other  by  comparing  the  RSTFs  and  analyzing  the  source  parameters  such  as 
source  duration,  source  dimension,  and  stress  drop.  We  conclude  that  the  EOF  method 
has  a  great  potential  for  discriminating  nuclear  explosions  and  earthquakes,  at  least  for 
events  with  mb  magnitudes  ranging  from  5.5  to  6.5,  by  comparing  their  RSTFs  and  source 

parameters. 
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Table  1.  Hypocentral  Parameters  of  Nuclear  Explosions  in  Central  Asia 


No. 

Date* 

Time,  UT 

Latitude, 

Longitude, 

Deptht, 

km 

mb 

Ms 

1 

881112 

0330:03.7 

50.078*  N 

79.988“  E 

0.30 

5.3 

Ig 

880614 

0227:06.4 

50.045*  N 

79.005“  E 

0.24 

5.0 

4.1 

2t 

880914 

0359:57.4 

49.833*  N 

78.808“  E 

0.50 

6.1 

4.5 

2g 

890708 

0346:57.6 

49.888“  N 

78.802“  E 

0.37 

5.6 

4.1 

3 

900816 

0459:57.6 

41.564*  N 

88.770“  E 

0.57 

6.2 

3g 

931005 

0159:56.5 

41.647*  N 

88.681“  E 

0.46 

5.9 

4.8 

4 

920521 

0459:57.5 

41.604“  N 

88.813“  E 

0.71 

6.5 

5.0 

4g-a 

931005 

0159:56.5 

41.647*  N 

88,681“  E 

0.46 

5.9 

4.8 

4g-b 

900816 

0459:57.6 

41.564“  N 

88.770“  E 

0.57 

6.2 

*Read  881112  as  Nov.  12,  1988:  year,  month,  day. 
fEstimated  with  the  formula  of  Jih  and  Wagner  [1991]. 
tJVE2. 
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Table  2.  Hypocentral  Parameters  of  Nuclear  Earthquakes  in  Central  Asia 


No. 

Date* 

Time,  UT 

Latitude, 

Longitude, 

Depth, 

km 

mb 

Ms 

1 

900515 

2229:59.3 

36.112"  N 

100.118"  E 

14 

5.5 

5.2 

Ig 

900428 

0356:46.7 

36.225"  N 

99.999"  E 

10 

4.8 

2 

900417 

0159:33.4 

39.436"  N 

74.900"  E 

33t 

6.0 

6.2 

2g 

900901 

2132:49.8 

39.375"  N 

74.760"  E 

17 

4.7 

3.8 

3 

900614 

1247:28.8 

47.869"  N 

85.076"  E 

58t 

6.1 

6.8 

3g 

900803 

0915:06.1 

47.963"  N 

84.961"  E 

33t 

6.0 

6.1 

4 

931002 

0842:32.8 

38.141"  N 

88.638"  E 

16 

6.2 

6.3 

4g 

931002 

0943:19.5 

38.127"  N 

88.502"  E 

12 

5.7 

5.3 

5 

930809 

1242:50.0 

36.348"  N 

70.840"  E 

233 

6.3 

5g 

930809 

1138:33.0 

36.395"  N 

70.707"  E 

230 

5.8 

6-1 

900426 

0937:10.9 

36.040"  N 

100.274"  E 

10 

5.7 

6-2 

900426 

0937:15.0 

35.986"  N 

100.245"  E 

8 

6.5 

6.9 

•  6-3 

900426 

0937:45.3 

36.239"  N 

100.254"  E 

10 

6.3 

6g 

900507 

0517:37.6 

36,032"  N 

100.341"  E 

33t 

5.3 

5.0 

7 

920819 

0204:37.4 

42.142"  N 

73.575"  E 

27 

6.6 

lA 

Vg 

920819 

0312:04.9 

42.107"  N 

73.264"  E 

21 

6.0 

6.3 

*Read  900515  as  May  15,  1990:  year,  month,  day. 
tDepths  are  not  well  constrained. 
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Table  3.  Source  Parameters  of  Nuclear  Explosions  in  Central  Asia 


No. 

Date* 

Ms 

10*^  Nm 

M.STF, 
10‘^  Nm 

'  /c 

Hz 

Re, 

km 

^CTt. 

MPa 

ig 

880614 

4.1 

5.0 

0.08 

1 

881112 

5.3 

0.17 

2.1 

0.17 

3.3±1.0 

0.25±0.10 

13 

(5-61) 

2g 

890708 

4.1 

5.6 

0.36 

2-1 

880914 

4.5 

6.1 

1.27 

3.4 

1.22 

2.5±0.6 

0.33±0.11 

41 

(18-140) 

2-2 

0.6 

0.22 

3g 

931005 

4.8 

5.9 

0.76 

3 

900816 

6.2 

1.63 

2.2 

1.67 

2.4±0.6 

0.34±0.12 

52 

(21-191) 

4g-a 

931005 

4.8 

5.9 

0.76 

4-1 

920521 

5.0 

6.5 

3.49 

4.3 

3.27 

1.6±0.4 

0.52±0.15 

28 

(13-  79) 

4-2 

1.2 

0.91 

4g-b 

900816 

6.2 

1.63 

4-1 

920521 

5.0 

6.5 

3.49 

1.9 

3.10 

1.6±0.4 

0.52±0.15 

27 

(13-  75) 

4-2 

0.6 

0.98 

*Read  880614  as  Junel4,  1988:  year,  month,  day. 

fValues  in  parentheses  are  the  minimum  and  maximum  value,  respectively. 
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Table  4.  Comparison  of  Spall  Slapdown  Phases  at  Different  Test  Sites 


Event 

Date* 

mfy 

Depth, 

km 

s 

Pss  Width, 

s 

Method/Reference/Remark 

Cannikin 

711106 

7.0 

1.79 

1.96 

0.6 

TD/Bakun  and  Johnson  [1973]/ATS 

1.95 

IC/Lay  et  al.  [1984al 

1.85 

SM/King  et  al.  [1974] 

0.7-1.4 

NF/[Springer  [1974]  and  Lay  [1991] 

0.8-1. 3 

LWP/Merritt  [1972] 

Milrow 

691002 

6.5 

1.22 

1.37 

0.5 

TD/Bakun  and  Johnson  [1973]/ ATS 

1.35 

IC/Lay  et  al.  [1984a] 

1.35 

SM/King  et  al.  [1974] 

Lop  Nor-92 

920521 

6.5 

0.7lt 

0.75 

0.8 

TEGF/this  study/XTS 

Boxcar 

680426 

6.3 

1.16 

1.35 

1.0 

NF/Springer  [1974J/NTS 

JVE2 

880914 

6.1 

0.50t 

0.45 

0.6 

TEGF/this  study/KTS 

Longshot 

651029 

6.0 

0.70 

1.40 

NF/Springer  [1974]/  ATS 

0.87 

SM/King  et  al.  [1974] 

Hardin 

870430 

5.8 

0.63 

1.30 

NF/Patton  [1990]/NTS 

Chancellor 

830901 

5.7 

0.62 

0.4-0.6 

NF/Patton  [1990]/NTS 

Baseball 

810115 

5.6 

0.56 

1.75 

NF/Taylor  and  Randall  [1989]/NTS 

Amarillo 

890627 

4.9 

0.60 

0.4-0.8 

NF/Stump  and  Reinke  [1991]/NTS 

Hardhat 

620215 

4.2 

0.29. 

1.40 

0.7 

NF/Patton  [1990]/NTS 

Borrego 

780213 

3.8 

0.56 

0.17 

NF/Taylor  and  Randall  [1989]/NTS 

Merlin 

650216 

0.30 

1.0-1. 2 

NF/Murphy  [1991]/NTS 

Hupmobile 

680118 

0.25 

0.45 

0.3 

NF/Springer  [1974]/NTS 

IC,  intercorrelation;  LWP,  lake  water  pressure;  NF,  near  field;  SM,  spectral  method; 
TD,  teleseismic  deconvolution;  TEGF,  teleseismic  EGF. 


ATS,  Amchitka  Test  Site;  XTS,  Xinjiang  Test  site;  NTS,  Nevada  Test  Site;  KTS,  Kazakhstan  Test  Site. 
*Read  711106  as  November  6,  1971:  year,  month,  day. 
fEstimated  with  the  formula  of  Jih  and  Wagner  [1991]. 
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Table  5.  Broadband  P  Waveform  Data  of  Earthquakes 
Used  in  This  Study 


No. 

Date*  Time,  UT 

Broadband  Seismic  Stations 

1 

Ig 

900515  2229:59.3 
900428  0356:46.7 

LZH 

2 

2g 

900417  0159:33.4 
900901  2132:49.8 

GAR  WMQ 

3 

3g 

900614  1247:28.8 
900803  0915:06.1 

BJI  KMI  LZH 

4 

4g 

931002  0842:32.8 
931002  0943:19.5 

AAK  ARU  KEV  KMI  LS  A  LZH 
TLYWMQ 

5 

5g 

930809  1242:50.0 
930809  1138:33.0 

ANTO  CHTO  COL  KONO 
MAJOPAB  TBT 

6 

6g 

900426  0937:15.0 
900507  0517:37.6 

KMI  WMQ 

7 

7g 

920819  0204:37.4 
920819  0312:04.9 

LZH  MDJ  WMQ 

* 

Read  900515  as  May  15, 

1990:  year,  month,  day 

Table  6.  Source  Parameters  of  Earthquakes  in  Central  Asia 


No. 

Date* * * § 

My  rtif, 

mCMT_ 

10‘'^Nm 

M„STF, 
10*’ Nm 

s 

r , 

km 

ZiCT, 

MPa 

Ig 

900428 

4.8 

1 

900515 

5.2  5.5 

1.1 

15.0 

1,5±0.1 

3, 5+0.2 

1.12 

(0,95-1. 34)t 
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Figure  1.  Map  of  Central  Asia  showing  the  study  area  and  epicentral  locations  of  nuclear 
explosions  (squares)  and  earthquakes  (circles).  Triangles  with  three-letter  codes  are  some 
of  the  broadband  seismic  stations  used  in  this  study.  The  numbers  marked  around  the 
earthquakes  correspond  to  those  listed  in  Tables  2,  5,  and  6. 


144 


3S11-2Q  KTS 

□ 

88061  4- 


890708 


(mMI) 

□  8809  14. 


0  km  5l 


(mM.9) 

931005 

□ 


□  9205211 


□ 

900816 


0  km  5l 


Figure  2.  (a)  PDE  locations  for  four  nuclear  explosions  at  the  Kazakhstan  Test  Site 
(KTS).  (b)  Relative  relocations  for  the  same  events  at  KTS  from  SPOT  satellite  image 
analysis  [Thurber  et  al,  1993],  (c)  PDE  locations  for  three  nuclear  explosions  at  the 
Xinjiang  Test  Site  (XTS).  (d)  Relocation  results  for  the  same  events  at  XTS  determined 
with  waveform  correlation  analysis  and  relative  location  technique. 
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Figure  3.  (a)  Vertical  component  broadband  P  waveforms  of  two  XTS  explosions  (mb  = 
6.5  and  6.2)  recorded  at  station  OBN,  Obninsk,  Russia,  36°  away  from  XTS.  The  P 
waveform  of  the  smaller  event  is  used  as  an  EOF.  The  numbers  below  the  traces  are 
negative  peak  amplitudes  in  digital  counts,  (b)  Retrieved  relative  STF  for  the  larger 
explosion.  Note  that  there  is  a  secondary  phase  on  the  seismogram  of  the  larger 
earthquake,  suggesting  the  secondary  pulse  of  RSTF  of  the  920521  explosion  is  indeed 
caused  by  the  explosion  source  rather  than  propagation  or  structure  effect. 
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Figure  4.  (a)  A  single-station  estimate  of  the  relative  STF  for  the  881112  explosion  (mb 
=  5.3)  in  KTS.  (b)  Single-station  estimates  and  a  stacked  trace  of  relative  STF  for  the 
880914  JVE2  explosion  (mb  =  6.1)  in  KTS.  Numbers  following  station  codes  are  azimuth 
and  epicentral  distance  in  degrees.  The  numbers  above  the  RSTF  traces  are  the  peak 
amplitudes  of  the  source  pulses. 
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Source  Time  Function:  900816  05:00  Loo  Nor  Explosion  mb=6.2 


Source  Time  Function:  920521  05:00  Loo  Nor  Explosion  mb=6.5 
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Figure  5.  (a)  The  single-station  estimated  and  stacked  RSTF  for  the  900816  explosion 
(mb  =  6.2)  in  XTS,  using  the  931005  explosion  (mb  =  5.9)  as  the  EGF  event,  (b)  Six 
single-station  estimates  of  RSTF  and  the  stacked  trace  for  the  920521  explosions  (mb  = 
6.5)  in  XTS,  using  the  931005  explosion  (mb  =  5.9)  as  the  EGF  event,  (c)  RSTF  of  the 
920521  explosion  in  XTS  estimated  with  the  900816  explosion  (mb  =  6.2)  as  the  EGF 
event.  Numbers  following  station  codes  are  azimuth  and  epicentral  distance  in  degrees. 
The  numbers  above  the  stacked  RSTF  traces  are  the  peak  amplitudes  of  the  source  pulses. 
Note  that  the  stacked  RSTFs  in  Figures  5b  and  5c,  determined  with  two  EGF  events,  are 
very  similar  except  for  the  peak  amplitudes. 
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Figure  6.  (a)  Convolution  of  an  assumed  STF  of  the  EGF  event,  Sg(f),  with  the  RSTF, 
Srit),  to  predict  an  "absolute"  STF,  S/(0,  for  the  920521  XTS  explosion  (m^  =  6.5).  The 
Sg(t)  are  assumed  to  be  triangular  pulses  with  pulse  widths  of  0.2,  0.4,  0.6,  0.8,  and  1.0 
s  and  the  peak  amplitudes  are  scaled  to  be  10,  5,  3.3,  2.5,  and  2  so  that  the  areas  under  the 
STFs  are  one  unit.  The  horizontal  bars  are  the  minimum  and  maximum  pulse  widths  of 
predicted  STFs.  (b)  Spectra  of  the  predicted  STFs  and  the  RSTF  for  the  920521  mb  =  6.5 
explosion.  Horizontal  bar  presents  variation  of  the  comer  frequency, /c,  of  the  RSTF  and 
the  predicted  STF. 
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Figure  7.  Comparison  of  RSTFs  for  the  four  explosions  at  KTS  and  XTS.  The 
numbers  in  parentheses  are  mb  magnitudes  for  the  EOF  events.  The  RSTF  of  two 
explosions  with  mb  =  5.3  and  6.2  show  a  simple  pulse.  The  other  two  explosions  (mb  = 
6. 1  and  6.5)  have  the  secondary  pulses  with  smaller  amplitudes,  which  we  explain  as  the 
spall  slapdown  phases. 
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Soum  Time  Function:  931002  08:42  Xinjiang  EartnouaKe  mD=6^ 


Source  Time  Function:  930809  12:42  Hindu  Kush  EarthouaKe  fnbs6.3 


Figure  8.  (a)  Eight  single  station  estimates  of  RSTF  for  the  931002  Xinjiang  earthquake 
(mb  =  6.2)  and  the  stacked  RSTF  trace,  indicating  a  small  precursor  event  before  the 
mainshock.  An  aftershock  {mb  =  5.7)  was  used  as  the  EGF  event,  (b)  Seven  single- 
station  estimates  of  the  RSTF  and  the  tacked  RSTF  trace  for  the  930809  Hindu  Kush  deep 
earthquake  {mb  =  6.3)  indicate  that  the  mainshock  is  a  double  event.  The  EGF  event  was  a 

=  5.8  foreshock. 
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Figure  9.  Comparison  of  RSTFs  for  the  seven  earthquakes  in  Central  Asia.  Top  trace  is 
a  single-station  estimate,  and  the  other  traces  are  stacked  RSTFs.  The  numbers  in  brackets 
are  mb  magnitudes  for  the  EGF  events.  The  locations  of  the  events  are  shown  in  Figure  1 
and  listed  in  Table  2.  Relative  source  durations  for  the  earthquakes  range  from  7  to  29  s. 
Most  of  RSTFs  show  the  earthquakes  have  a  complex  rupture  process.  The  RSTF  of  the 
900426  Qinghai  mainshock  sequence  indicated  that  it  started  with  a  smaller  precursor  event 
about  5  s  before  the  second  subevent,  followed  by  the  third  subevent  about  30  s  later.  PDE 
reported  the  origin  time  of  the  three  subevents  was  at  0937:10.9,  0937:15.0,  and 
0937:45.3  (Table  2). 
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Figure  10.  (a)  Convolution  of  an  assumed  STF  of  the  EOF  event,  Sg(t),  with  the  RSTF, 
Sr(t),  to  retrieve  an  "absolute"  STF,  5/(0,  for  the  Kazakh  mainshock  (mb  =  6.1).  The 
Sg(t)  are  assumed  to  be  triangular  pulses  with  source  duration  of  9  ±  1  s,  and  the  peak 
amplitudes  are  scaled  to  be  1/4,  9/2,  and  1/5  so  that  the  areas  under  the  STFs  equal  one. 
(b)  Estimate  of  the  rise  time  from  the  three  "absolute"  STFs  of  the  larger  Kazakh  event,  (c) 
Recovery  of  an  "absolute"  STF  for  the  Hindu  Kush  mainshock.  The  source  pulse  widths 
of  the  EOF  event  are  assumed  to  be  6,  7,  and  8  s,  and  the  STF  peak  amplitudes  are  scaled 
to  be  1/3,  2/7,  and  1/4,  respectively,  (d)  Estimate  of  the  rise  time  for  the  first  subevent 
from  the  three  "absolute"  STFs  of  the  Hindu  Kush  mainshock.  The  horizontal  bars  in 
Figures  10b  and  lOd  represent  errors  of  the  rise  time  measurements. 
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Figure  11.  Comparison  of  relative  source  durations  for  earthquakes  and  explosions  in 
Central  Asia.  Note  that  the  source  durations  for  explosions  are  typically  less  than  those  of 
earthquakes  by  a  factor  of  10. 
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ABSTRACT 

The  discrimination  of  seismic  events  as  small  as  magnitude  2.0  at  regional  distances  is  of  critical 
importance  to  achieving  the  goal  of  monitoring  the  Comprehensive  Test  Ban  Treaty  (CTBT). 
Although  many  discriminants  work  well  for  large  earthquakes  and  explosions,  little  has  been  learned 
about  their  performance  for  smaller  events.  In  this  study,  we  explore  the  performance  of  Pg/Lg 
and  Lg  spectral  ratio  discriminants  on  smaller  events  using  waveforms  of  14  earthquakes  and  10 
industrial  explosions  with  magnitudes  between  1.6  and  3.8  in  New  England,  recorded  by  short- 
period  stations  at  distances  up  to  250  km.  We  find  that  the  Pg/Lg  discriminant  fails  to  separate  the 
small  explosion  and  earthquake  populations,  but  Lg  spectral  ratio  appears  to  be  a  good  discriminant 
for  these  events.  Our  analysis  indicates  that  the  Lg. spectra  observed  from  industrial  explosions 
have  much  larger  amplitudes  at  lower  frequencies,  while  earthquakes  have  higher  amplitudes  at 
higher  frequencies,  which  is  reversed  from  the  situation  for  large  earthquakes  and  explosions.  To 
understand  the  physical  basis  of  the  Lg  spectral  ratio  discriminant  for  small  magnitude  seismic 
events,  we  performed  numerical  modeling  experiments  to  study  the  effects  of  source  properties, 
depth,  and  attenuation.  Synthetic  seismograms  demonstrate  that  the  lower  frequency  content  of 
small  explosions,  relative  to  earthquakes,  is  partially  explained  by  their  shallower  depth.  The 
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energy  from  a  near-surface  explosion  experiences  greater  high  frequency  attenuation  in  shallow, 
low  Q,  layers  and  greater  trapping  of  low  frequency  energy  in  lower  velocity  layers.  Furthermore, 
the  separation  between  earthquake  and  explosion  Lg  spectra  is  enhanced  by  two  source  effects: 
the  shorter  source  duration  of  small  earthquakes  and  the  extended  duration  of  explosions  that  are 
multiply  detonated.  These  modeling  results,  together  with  our  empirical  study  of  New  England 
events,  suggest  that  Lg  spectral  ratio  is  a  promising  method  for  discriminating  seismic  events  of 
small  magnitude. 


OBJECTIVES 

The  objective  of  this  project  is  to  characterize  the  source  properties  of  seismic  events,  as  a  function 
of  event  type  and  magnitude,  with  the  ultimate  goal  of  developing  effective  regional  discriminants. 
In  previous  work  we  applied  the  EGF  method  to  estimate  source  time  functions  for  a  wide  variety  of 
events  (earthquakes,  chemical  and  nuclear  explosions,  mining  tremors)  over  a  large  magnitude  range 
(0. 9-6.6),  and  discovered  that  the  source  time  function  duration  could  discriminate  earthquakes  and 
explosions  in  certain  magnitude  ranges,  including  very  small  events  (magnitudes  of  3.0  or  less)  but 
not  including  events  in  the  magnitude  range  3-4.5. 

To  achieve  the  goal  of  monitoring  a  CTBT,  it  is  critical  to  discriminate  seismic  events  as  small 
as  magnitude  2.0  (Murphy,  1995).  Although  many  methods  work  well  for  discriminating  larger 
earthquakes  and  explosions,  little  has  been  learned  about  their  performance  on  the  smaller  events. 
The  objective  of  our  present  work  is  to  understand  the  physical  basis  of  our  previous  results  for 
small  events,  and  to  compare  them  to  the  behavior  of  discriminants  based  on  Lg. 

RESEARCH  ACCOMPLISHED 

In  our  previous  studio  (Toksoz  et  al.,  1993;  Li  et  al.,  1994,  1995a, b)  we  applied  an  empirical 
Green’s  function  (EGF)  method  to  characterize  seismic  sources  and  estimate  source  time  functions 
for  explosions,  earthquakes,  and  mining  tremors  covering  the  magnitude  range  0.9-6.6.  Figure  1 
displays  source  duration  vs.  event  magnitude  for  all  the  events  we  have  studied  .to  date.  We  see 
that  the  earthquake  and  explosion  populations  appear  to  obey  different  source  scaling  relations. 
For  magnitudes  between  3.5  and  4.5,  the  earthquake  and  explosion  populations  overlap  due  to  the 
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crossover  of  two  scaling  relations  with  different  slopes,  indicating  that  simple  source  duration  is  not 
a  good  discriminant  in  this  magnitude  range.  However,  for  M  =  5  and  greater,  source  durations  of 
nuclear  explosions  are  about  a  factor  of  5  to  10  shorter  thatt  those  of  the  earthquakes  (Figure  1). 
Thus  source  duration,  like  other  discriminants  such  as  P/Lg,  Lg  spectral  ratio,  and  mb-Ms  (e.g., 
Pomeroy  et  al.,  1982;  Bennett  and  Murphy,  1986;  Taylor  et  al,  1988),  works  as  a  good  discriminant 
for  larger  earthquakes  and  explosions. 

For  the  small  (M<3)  events  we  analyzed,  the  discriminant  also  works  but  in  an  opposite  way. 
The  source  durations  of  quarry  and  construction  blasts  are  significantly  longer  than  those  of  earth¬ 
quakes  with  similar  magnitudes  (Figure  1).  This  is  consistent  with  the  results  of  Herrin  et  al.  (1994), 
who  found  that  Lg  spectra  of  small  magnitude  earthquakes  in  the  Vogtland  area  have  relatively 
higher  peak  frequencies  compared  to  explosions  with  similar  magnitudes.  Possible  explanations 
for  these  observations  from  small  earthquakes  and  explosions  are  source,  depth,  and  attenuation 
effects,  i.e., 

1.  Smaller  earthquakes  and  explosions  may  have  significantly  different  source  dimensions  over 
space  and  time.  For  example,  the  source  duration  of  industrial  blasts  is  often  controlled  by 
shooting  practices,  such  as  multiple  detonation  and  ripple  firing. 

2.  Explosions  typically  are  detonated  at  shallow  depths  (<1  km)  while  earthquakes  occur  at 
depths  of  a  few  kilometers  or  more. 

3.  Explosions  are  typically  in  a  near-surface  medium  with  high  attenuation  while  earthquakes 
occur  in  a  medium  with  lower  attenuation. 

In  this  study  we  explore  the  performance  of  the  P/Lg  discriminant  and  Lg  spectral  ratio  dis¬ 
criminant  on  the  smallest  events  from  our  previous  work.  These  include  earthquakes  and  industrial 
explosions  with  magnitudes  between  1.6  and  3.8  in  Nev/  England  recorded  by  short-period  stations 
of  the  M.I.T.  seismic  network  at  regional  distances. 

Figure  3  shows  velocity  seismograms  of  one  quarry  blast  (M=2.0)  and  one  earthquake  (M=2.8) 
at  Littleton,  Massachusetts,  recorded  at  station  DNH,  which  is  about  80  km  from  the  epicenters 
(see  Figure  4).  As  we  expected,  the  small  earthquake  has  a  much  richer  high  frequency  content 
compared  to  the  explosion.  Bandpass  filtered  seismograms,  using  2-4  and  6-8  Hz  bands,  are  also 
shown  in  Figure  3.  For  the  explosion,  the  Lg  peak  amplitude  on  the  2-4  Hz  trace  is  about  5 
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times  greater  than  that  on  the  6-8  Hz  seismogram.  In  contrast,  the  Lg  peak  amplitude  on  the 
2-4  Hz  seismogram  is  about  2  times  smaller  than  that  on  the  6-8  Hz  seismogram.  The  Fourier 
spectra  of  the  Lg  phases  from  the  filtered  traces  (bottom  of  Figure  3)  show  this  same  pattern, 
with  the  explosion  displaying  relatively  more  low  frequency  energy.  Since  the  two  events  have  very 
similar  epicenters  and  are  recorded  at  a  common  station,  we  believe  that  the  observed  differences 
between  their  Lg  spectral  patterns  are  mainly  caused  by  differences  between  the  sources  themselves 
and  by  effects  related  to  focal  depth  differences,  as  opposed  to  propagation  along  different  paths. 
This  suggests  that  Lg  spectral  ratios  may  generally  discriminate  between  small  earthquakes  from 
explosions. 

It  is  interesting  to  compare  this  result  to  an  example  with  large  events.  Figure  2  shows  Lg 
spectra  of  one  Xinjang  earthquake  (931002,  77156.2)  and  one  Lop  Nor  explosion  (900816,  77156.2) 
recorded  at  station  ARU  about  25  degrees  away.  We  see  that  the  Lg  spectral  ratio  between  higher 
(0.4  to  2  Hz)  and  lower  (0.04  to  0.2  Hz)  frequency  bands  is  much  higher  for  the  explosion  than 
for  the  earthquake.  That  is,  the  large  explosion  has  relatively  higher  peak  frequencies  compared 
to  the  large  earthquake,  which  is  reversed  from  the  small  event  example  of  Figure  3.  Again,  the 
two  events  have  similar  propagation  paths  and,  considering  also  the  lower  frequencies  used,  the 
Lg  spectra  are  probably  dominated  by  the  source  effect.  The  behavior  of  Lg  spectral  ratios  for 
small  vs.  large  events  suggests  that  the  source  scaling  relations  for  source  duration  may  be  true  in 
general. 

Next,  we  apply  the  Lg  spectral  ratio  discriminant  to  10  explosions  (M=1.6  to  2.7)  and  14 
earthquakes  (M=2.2  to  3.8)  in  New  England  (see  Figure  4).  The  spectral  ratio  Lg(2-4Hz)/Lg(6- 
8Hz)  plotted  against  epicentral  distance  and  magnitude  are  shown  in  Figure  5.  We  find  that  the 
Lg  spectral  ratio  separates  well  the  earthquake  and  explosion  populations,  although  the  magnitude 
and  distance  ranges  containing  both  types  of  events  are  limited.  We  note  that  Pg/Lg  at  either  2-4 
Hz.  or  6-8  Hz  does  not  separate  the  small  explosion  and  earthquake  populations. 

To  understand  the  physical  basis  for  the  Lg  spectral  ratio  discriminant  for  small  magnitude 
events,  we  performed  tests  with  synthetic  data  to  study  the  various  effects  listed  above:  source 
properties,  depth,  and  attenuation.  Synthetic  seismograms  were  computed  using  the  discrete 
wavenumber  algorithm  of  Mandal  and  Toksoz  (1990)  and  a  1-D  structure  comprising  three  crustal 
layers  over  a  half-space  (Table  1).  Seismograms  were  computed  for  an  explosion  source  at  0.1  km 


158 


depth  and  a  thrust  earthquake  at  a  depth  of  4  km.  The  source-receiver  distance  is  80  km.  We 
selected  both  source  time  functions  to  be  a  triangle  pulse  with  time  duration  of  0.1  s.  The  synthetic 
velocity  seismograms  for  the  explosion  and  earthquake  are  shown  in  Figures  6b  and  6c,  and  the  Lg 
spectra  are  displayed  in  Figures  7b  and  7c.  From  the  seismograms  and  Lg  spectra  we  can  see  that, 
for  the  same  source  pulse,  the  explosion  has  much  higher  amplitudes  at  lower  frequencies  (<  4Hz) 
while  the  earthquake  has  higher  amplitudes  at  higher  frequencies  (>  6Hz).  Since  the  explosion  is 
detonated  in  the  more  attenuating  surface  layer,  and  ray  paths  go  through  this  layer  twice,  there 
is  a  greater  loss  of  high  frequency  energy  from  this  source  compared  to  the  earthquake,  which  is  in 
the  second  layer  of  the  model.  There  may  also  be  an  enhancement  of  lower  frequency  energy  from 
the  explosion  due  to  the  trapping  of  energy  in  the  shallow  layer. 

We  also  calculated  synthetics  for  a  multiply  detonated  blast  at  a  depth  of  0.1  km  (Figures  6a  and 
7a).  The  multiple  explosion  has  10  subevents  with  delay  times  varying  between  0.08  to  0.12  s  and  a 
total  source  duration  of  1.0  s.  We  see  that  multiple  firing  increases  the  energy  at  lower  frequencies 
and  decreases  it  at  higher  frequencies,  thus  increasing  the  difference  from  the  earthquake  synthetic. 

As  shown  in  Figure  1,  the  source  duration  for  a  small  earthquake  is  typically  2  to  3  times 
shorter  than  for  an  explosion  with  similar  magnitude.  To  examine  how  a  sharper  source  pulse  for 
earthquakes  affects  the  Lg  ratio  discriminant,  we  recalculated  the  earthquake  synthetic  with  the 
triangle  source  pulse  reduced  to  0.05  s  in  duration  (Figures  6d  and  7d).  As  expected,  the  energy 
content  shifts  to  higher  frequency,  further  enhancing  the  difference  between  the  earthquake  and 
explosions. 

We  applied  the  Lg(2-4Hz)/Lg(6-8Hz)  discriminant  to  the  synthetic  data  from  all  four  cases. 
The  Lg  spectra  in  the  two  frequency  bands  are  shown  for  the  multiple  blast,  single  blast,  and 
earthquakes  with  differing  source  durations  in  Figures  8a-8d,  respectively.  The  Lg  spectral  ratio 
between  the  two  frequency  bands  clearly  separates  the  explosions  and  earthquakes.  The  greatest 
separation  is  between  the 'multiple  explosion  (Figure  8a)  and  the  earthquake  with  a  shorter  source 
duration  (Figure  8d). 
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CONCLUSIONS  AND  RECOMMENDATIONS 


The  Lg  spectral  ratio,  using  frequency  bands  of  2-4  Hz  and  6-8  Hz,  discriminates  earthquakes  and 
industrial  explosions  in  New  England  at  magnitudes  below  3  and  epicentral  distances  up  to  250  km. 
The  Pg/Lg  discriminant  is  not  effective  on  these  events,  however.  The  success  of  the  Lg  spectral 
ratio  discriminant  is  consistent  with  our  previous  work  that  showed  that  the  source  time  duration 
of  M<3  earthquakes  is  2  to  3  times  shorter  than  that  of  explosions  of  similar  magnitude. 

Based  on  numerical  modeling  of  earthquake  and  explosion  seismograms  in  New  England,  we 
conclude  that  the  lower  frequency  content  of  small  explosions,  relative  to  earthquakes,  is  partially 
attributable  to  their  shallower  depth.  The  energy  from  near-surface  explosions  experiences  greater 
high  frequency  attenuation  in  shallow,  low  Q  layers  and  greater  trapping  of  low  frequency  energy 
in  lower  velocity  layers.  The  separation  between  earthquake  and  explosion  Lg  spectra  is  enhanced 
by  two  source  effects:  the  shorter  source  duration  of  small  earthquakes  and  the  extended  duration 
of  explosions  that  are  multiply  detonated. 

These  results  imply  that  source  time  duration  and  Lg  spectral  ratio  are  promising  discriminants 
for  small  earthquake  and  explosions.  The  Lg  discriminant  is  particularly  interesting  since  Lg  spectra 
are  more  directly  and  easily  measured.  For  future  work,  therefore,  we  recommend  that  Lg  spectral 
discriminants  be  tested  on  a  large  number  of  events  from  different  regions  to  assess  their  general 
effectiveness  and  transportability.  Of  particular  importance  is  to  learn  their  effectiveness  on  events 
of  magnitude  3  and  greater,  and  at  greater  distances.  To  better  understand  the  physical  basis  for 
the  Lg  spectral  ratio  discriminant,  and  to  provide  the  knowledge  needed  to  adapt  the  discriminant 
to  different  regions,  we  recommend  further  numerical  modeling  studies  aimed  at  determining  how 
Lg  spectra  are  affected  by  vertical  layering  in  the  crust,  and  2-D  and  3-D  scattering  by  topographic 
and  geologic  features. 


160 


REFERENCES 


Bennett,  T.,  and  J.  Murphy,  1986.  Analysis  of  seismic  discrimination  using  regional  data  from 
western  United  States  evnets,  Bull.  Seis.  Soc.  Am.,  76,  1069-1086. 

Herrin,  E.,  Burlacu,  V.,  Gray,  H.L.,  Swanson,  J.,  Golden,  P.,  and  Myers,  B.,  1994.  Research  in 
regional  event  discrimination  using  Ms:mb  and  autoregressive  modeling  of  Lg  waves.  Proceed¬ 
ings,  16th  Annual  Seismic  Research  Symposium,  Thornwood,  New  York,  Phillips  Laboratory, 
Hanscom  AFB,  Massachusetts,  152-158.  PL-TR-94-2217. 

Li,  Y.,  Rodi,  W.,  and  Toksoz,  M.N.,  1994.  Seismic  source  characterization  with  empirical  Green’s 
function  and  relative  location  techniques,  Proceedings,  16th  Annual  Seismic  Research  Sym¬ 
posium,  Thornwood,  New  York,  Phillips  Laboratory,  Hanscom  AFB,  Massachusetts,  231-237. 
PL-TR-94-2217 

Li,  Y.,  Toksoz,  M.N.,  and  Rodi,  W.,  1995a.  Source  time  functions  of  nuclear  explosions  and  earth¬ 
quakes  in  central  Asia  determined  using  empirical  Green’s  functions,  J.  Geopbys.  Res.,  100, 
659-674. 

Li,  Y.,  Rodi,  W.,  and  Toksoz,  M.N.,  1995b.  Discrimination  of  earthquakes,  explosions,  and  min¬ 
ing  tremors  using  the  empirical  Green’s  function  method.  Proceedings,  17th  Annual  Seismic 
Research  Symposium,  Scottdale,  Arizona,  Phillips  Laboratory,  Hanscom  AFB,  Massachusetts, 
78-87.  PL-TR-95-2108 

Mandal,  B.  and  Toksoz,  M.N.,  1990.  Computation  of  complete  waveforms  in  general  anisotropic 
media  -  results  from  an  explosion  source  in  an  anisotropic  medium,  Geopbys.  J.  Int.,  103,  33-45. 

Murphy,  J.  R.,  1995.  An  overview  of  seismic  Discrimination  issues  relevant  to  CTBT  monitoring. 
Proceedings,  17th  Annual  Seismic  Research  Symposium,  Scottdale,  Arizona,  Phillips  Labora¬ 
tory,  Hanscom  AFB,  Massachusetts,  88.  PL-TR-95-2108 

Pomeroy,  P.W.,  J.B.  Best,  and  T.V.,  McEvilly,  1982.  Test  ban  treaty  verification  with  regional  data 
-  a  review.  Bull  Seis.  Soc.  Am.,  72,  S^9-S129. 

Talyor,  S.R.,  N.W.Sherman,  and  D.D.  Marvin,  1988.  Spectraal  discrimination  between  NTS  ex¬ 
plosions  and  western  United  States  earthquakes  at  regional  distances.  Bull.  Seis.  Soc.  Am.,  78, 
1563-1579. 

Toksoz,  M.N.,  Li,  Y.,  and  Rodi,  W.,  1993.  Seismic  source  characterization  with  empirical  Green’s 
function  and  relative  location  techniques.  Proceedings,  15th  Annual  Seismic  Research  Sympo- 


161 


Magnitude 


Figure  1:  Source  duration  vs.  magnitude  for  earthquakes  (circles),  explosions  (squares),  and  mining 
tremors  (triangles),  determined  by  EGF  analysis  (Li  et  al.,  1995b).  For  events  having  magnitude 
less  than  3  or  greater  than  4.5,  the  source  duration  is  a  good  discriminator  between  earthquakes  and 
explosions.  In  the  remaining  magnitude  range  (3.0  to  4.5)  the  earthquake,  explosion  and  mining 
tremor  populations  overlap. 

Xinjiang  Earthquake  Vs.  Lop  Nor  Nuclear  Explosion 


Figure  2:  The  Lg  spectra  for  an  77156.2  earthquake  in  Xinjang  (931002)  and  an  77156.2  Lop  Nor 
nuclear  explosion  (900816)  recorded  at  broadband  station  ARU. 
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Figure  3:  Short  period,  vertical  component  seismograms  recorded  at  station  DNH  (distance  80  km) 
for  a  Littleton  explosion  (top)  and  a  nearby  earthquake  (middle).  Unfiltered,  bandpass  filtered  from 
2-4  Hz,  and  bandpass  filtered  from  6-8  Hz  waveforms  are  displayed  in  velocity.  Lg  spectra  of  the 
two  frequency  bands  are  plotted  for  the  explosion  (bottom  left)  and  the  earthquake  (bottom  right). 
Note  that  the  Lg  spectral  ratios  between  the  two  frequency  bands  are  different  for  the  explosion 
and  the  earthquake. 
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Figure  4;  Map  of  New  England  showing  locations  of  10  explosions  (crosses)  and  14  earthquakes 
(circles)  and  5  M.I.T.  seismic  stations  (triangles)  used  in  this  study.  The  distance  from  the  Littleton 
earthquakes  and  explosions  (close  to  station  WFM)  to  station  DNH  is  about  80  km. 

Lg  Spectral  Ratio  Vs.  Distance  Lg  Spectral  Ratio  vs.  Magnitude 


Distanced™)  Magnitude 


Figure  5:  The  Lg  spectral  ratio  (2-4Hz)/(6-8Hz)  versus  epicentral  distance  (left)  and  magnitude 
(right)  for  explosions  (crosses)  and  earthquakes  (circles).  The  Lg  spectral  ratio  well  separates  the 
small  earthquake  and  explosion  populations. 
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Figure  6:  Synthetic  velocity  waveforms  for  a  multiple  explosion  (a),  a  single  explosion  (b),  and  a 
thrust  earthquakes  with  a  source  duration  of  0.1  s  (c)  and  0.05  s  (d).  The  left  column  shows  the 
source  time  functions  used  in  the  calculation  of  the  synthetics. 


Frequency  (Hz) 


Figure  7:  Fourier  spectra  of  the  Lg  phase  for  a  multiple  explosion  (a),  a  single  explosion  (b),  and 
two  thrust  earthquakes  (c  and  d)  with  different  source  durations..  Note  that  the  explosions  are 
richer  in  lower  frequency  energy  and  earthquakes  have  more  higher  frequency  energy. 
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